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The following Notation has been used in this Thesis 
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SYNOPSIS 


The dynamic soil-structure interaction deals with response of structures 
and/or ground subjected to time-dependent loads. The loads could be 
applied directly on to the structure (arising from rotating machinery, 
wind, etc ) or thiough soil in the form of incident waves (arising out 
of earthquake or an explosion or unbalanced machmeiy etc. ). 
Analytical methods available for dynamic soil structure inteiaction arc 
applicable only to highly idealized geometric representations and to 
linear elastic material behaviour in frequency domain. Numerical 
methods which need domain discretization can be not be applied 
directly as soil exists as a half-space. The discretization of finite 
(computational) domain coupled with boundary conditions on the 
finite domain to account for energy radiation (approximately) can be 
used to model the problems. 

The Present investigation addresses the problems of dynamic soil- 
structure interaction in time domain. The computational schemes 
involves the use of FEM to discretize a large computational domain 
and a load dependent Lanezos transformation procedure to reduce the 
number of degrees of freedom obtained from the FEM. The resulting 
dynamic equilibrium equations are integrated using an adaptive strat- 
egy devised for a second order integration method. The computational 
scheme developed has the following capabilities — 1. No need for 
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energy absorbing boundaries/ 2 Retains all the advantages of MiM, 
3. An integrated analysis of soil, structure and foundation is possi- 
ble, 4 Can analyze all modes simultaneously, 5. Can handle several 
arbitrary nature, 6 Can be used for transient as well as for steady 
state responses, 7. Error estimates and adaptive error control strategies 
are available for temporal discretization as well as for transformation 
process. 

The study covers the areas of flexible and rigid foundation i espouse 
to various types excitations. Effect of shape and embeddment of the 
foundation have been investigated Soil has been been idealized as 
linear isotropic continuum. Non-homogenity of the soil has been 
considered in terms of continuously varying modulus and as a layered 
system Response of foundation on finite thick soil stratum has also 
been investigated. 

Response of pile foundation soil system has been investigated for vari- 
ous modes of excitation. Effect of various parameters like embeddment 
of pile cap, slenderness ratio relative rigidity and spacing of Individual 
piles have been studied 

Vibration isolation has been investigated with open and filled trenches. 
Various parameters affecting the amplitude reduction like material 
properties of the trench fill, the geometry and location of the trench 
and excitation frequency. 

In Chapter 1 a brief introduction to the problems is presented. The 
literature available is discussed and the need and scope of the present 
investigation is highlighted, Chapter 2 introduces the computational 
model utilized in this investigation, The validity of the elastic con- 
tinuum theory of the soil has been discussed. A brief inlrcduclion 
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to finite element method has been presented, and the application of 
the finite element method to dynamic problems is presented A laige 
computational domain has been used to represent semi-infinite soil 
medium, and the system has been discietized. Load dependent Lanc- 
zos vectors have been utilized for coordinate transformation Apriori 
error estimates have been used to asses the quality of the coordinate 
transformation and number of vectors or ordei of reduced system have 
been estimated by specifying error tolerance. The resulting dynamic 
equilibrium equations in reduced space have been solved in time do- 
main. A valiant of Newmark 0 method called 0 method, which lelams 
second order accuracy with controllable numerical damping has been 
utihzed An adaptive time stepping strategy has been devised by post 
processing the solutions with specihed local truncation error tolerance. 
The solutions (responses of required degrees of freedom ) are obtained 
in original space. 

An extensive validation has been carried out to verify the accuracy 
of the numencal model in hme domain. An half-space and half- 
space with vertical trench have been chosen. To simulate transient 
excitation Ricker's Wavelet loading has been prescribed. Both vertical 
and horizontal excitations have been investigated in two dimensions. 
Excellent agreement with published results is observed. Effect of 
number of Lanezos vectors, model size, number of elements on the 
accuracy of the results has been brought out. Different spatial variation 
of the loads have been considered for generation of Lanezos vectors 
and found that reduced system matrices can be used to analyze the 
problems whose spatial loading patterns are similar to the one used 
to generate Lanezos vectors. Computational times required for each 
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stage of the computation is reported for different number of Lanczos 
vectors. The present model has been found to represent non convex 
domains in time domain with good accuracy The solutions obtained 
are found to obey causality conditions for half-space and half-space 
with vertical trench problems. 

In Chapter 3, results of the detailed parametric studies carried out on 
embedded foundations have been presented in terms of the compliance 
functions. The results obtained with present model are compared with 
existing time-domain boundary element model solution for vertical, 
horizontal, rotational and torsional modes The agreement is found 
to be sahsfactory. Effect of depth of embeddment, aspect ratio of the 
footing, relative flexibility on the compliance functions for different 
modes at different frequency ratios have been investigated. Behaviour 
of flexible foundahons at different embeddmenls are presented to bring 
out the change in the response. Effect of linearly varying modulus, 
presence of rigid bedrock at finite depth and presence of layeiing in 
the soil system have been investigated. It has been found that effect 
of these parameters differ from verheal mode to horizontal mode of 
vibration. Linearly increasing modulus shows stronger effect on hoii- 
zontal mode than vertical mode. It increases the resonant fiequency 
for horizontal mode of vibration. Presence of hard rock is strongly 
influences the vertical compliance functions. The compliance functions 
are influenced by small thickness (usually 10 times the halfwidth of 
the foundation) of the top layer of the soil. As frequency of excita- 
tion increased the thickness of the soil layer influencing the behaviour 
decreases. 

In Chapter 4 results obtained for pile foundations in homogeneous and 
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layered soils have been presented Analysis of piles and pile cap has 
been earned out in an integrated manner along with soil medium The 
study has been carried out in such a way that the effect of presence 
of piles in foundation system has been brought out The responses 
in different modes of vibration in terms of compliance functions for 
pile cap are presented A few parametric studies on the behaviour of 
single piles are carried out to. The response of single piles foi transient 
excitation is investigated for validation and the behaviour predicted 
by the present model agrees well with those repoited. 

In case of single piles, increase in slenderness ratio substantially reduces 
compliance functions for vertical mode of vibration Whereas for 
horizontal and rotational modes of vibration the effect of slenderness 
ratios is felt at higher modular ratios only. 

The effect of group of piles on compliance functions of pile foundation 
(of the pile cap) is different in different modes of vibration. The pres- 
ence of piles (four) reduces the compliance functions for vertical and 
rotational modes substantially (by a factor of 4), whereas the reduction 
is of the order of 1.5 for horizontal mode of vibration, and by a factor 
of 10 for rotational modes. The effect of modular ratio on compli- 
ance functions of the pile cap has been investigated. It is found that 
modular ratio does not a have significant effect on vertical and hori- 
zontal modes, but has strong influence on rotational response. Effect 
of stiffer stratum below has been investigated. The presence of stiffer 
soil influences the vertical mode strongly and increases the resonant 
frequencies. Unlike the block foundation considerable thickness of top 
layer influences the pile group behaviour. Increase in number of piles 
(reducing the pile spacing) reduces the compliances considerably in 
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verhcal and rotahonal modes. 

A parametric study on machine pile foundation soil system (MPFS) 
transmitting periodic but not harmonic excitations simultaneously is 
presented to study real life situations. The parameters considered 
for the study are the slenderness ratio of the piles, shear modulus 
of the soil and pile cap/block thickness. The machine consists of a 
diesel generating set, and is transmitting horizontal force and rocking 
moment with one period, and torsional moment with another period. 
These excitations are arbitrary except that they are periodic individually 
The maximum displacements at the center of gravity of the foundation 
has been used to present the results. 

Increasing the shear modulus reduces the response drastically at low 
shear moduli!. The responses do not significantly reduce due to 
increase of the modulus beyond a value of QMPa. Presence of piles 
reduces the horizontal displacement by about 60% at low modulii and 
by a small^ fraction at higher modulii. For rotational and torsional 
modes the reduction is of the order of 100% due to presence of piles 

t 

at low modulii. Increasing the slenderness ratio beyond a value of 
20 in horizontal and rotational modes, and beyond a value of 10 in 
torsional mode, does not reduce the responses. The pile cap/block 
thickness has very little effect on horizontal and rotational modes. 
Increasing' the pile cap thickness beyond 1.0m increases the responses 
in some of the situations. 

In Chapter 5 a brief study has been carried out on the vibration 
isolahon systems. Open and filled trenches as well as concrete wall 
barriers are considered for the study. The effect of width, depth and 
location of the trench/barrier on vibration isolation efficiency, in teims 
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of amplitude reduction, has been carried out for vertical and horizontal 
modes Effect of layenng of soil has also been investigated. 

In Chapter 6 conclusions drawn from the present study aie summa- 
rized The highlights of the computational scheme viz the advantages 
of rahonal and integrated analysis for tiansient as well as steady state 
responses with small computational resources are reported. Possible 
extensions of the present scheme for non-linear analysis have been 
indicated 
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INTRODUCTION 


The dynamic soil-structure interaction deals with response of sti uctures 
and/or ground subjected to a specified time-varying loads The loads 
could be applied directly on to the structure (arising from rotating 
machinery, wind etc.) or through soil in the form of incident wave 
(arising out earthquake or an explosion). 

Dynamic soil-structure interaction includes analysis of machine foun- 
dation soil systems, which is an important area of study. The Dynamic 
analysis of foundation soil systems have received a lot attention of 
researchers recently. This is partly attributed to recognition of the 
importance of earthquake hazards which warrants dynamic analysis 
of structural systems including substructures (foundations), and avail- 
ability of digital computers with lot of memory and computing power 
required for dynamic analysis. The dynamic soil-structure interac- 
tion unlike structural dynamics counterpart is quite complicated and 
involved because of its infinite extent (domain) and uncontrollable 
variations of the characteristics (non homogeneity etc. ). Any realistic 
analysis of dynamic soil structure interaction needs three dimensional 
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idealization of the geometry. 

Analytical methods available for dynamic soil structure inleraclion aie 
applicable only to highly idealized geometric representations and to 
linear elastic material behaviour in frequency domain. Numerical 
methods which need domain discretization can not be applied diioctly 
as soil exists as a half-space The discretization of finite (computa- 
tional)domain coupled with appropriate boundary conditions on the 
hnite domain to account for energy radiation can be used to model 
the problems. 

The study of dynamic soil-structure interaction can be classilied in 
to two distinct approaches — frequency domain approach and lime 
domain approach. The frequency domain approach is suitable for 
linear elastic problems with periodic excitations. The lime domain 
approach can handle all type of excitations and non-linear material 
behaviour It also provides responses in a natural way ( i.c. as events 
occur). 

1.1 Brief Review of Existing Literature 

Early investigations of soil-structure interaction deall with re.spon.se 
of soil idealizing it as linear, homogeneous, semi-infinite ela.slic half- 
space. Lamb (1904) initiated study of response of elastic hnlf-space 
subjected to dynamic loads. Reissner (1936) making use of Lamb'.s 
solution, obtained the response of harmonic disc load on Ihu .surface 
of an elastic half-space. He assumed a uniform dislribulion of .slre.sse.s 
under the disc. Quinlan (1953) derived the equation.^ for obtaining 
the responses with contact pressures which vary across llie diameter 
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with parabolic distribution, with uniform distnbution and with the 
distiibution corresponding to a rigid base. He gave solutions for rigid 
base approximations Sung (1953) developed solutions for ciiculai 
contact area with three pressure distributions, namely, 1, uniform, 2 
parabolic, and 3. pressuie distribution corresponding to rigid base 
He gave solutions for different Poisson's ratios These solutions can 
be used to compute the compliance functions of rectangular footing 
making use of the equivalent circular area. This equivalent area method 
provides solutions of reasonable accuracy for aspect ratios up to two 
(P E Richart, 1960) He also developed mathematical expressions for 
rectangular footing resting on half-space. Kobori (1962), Thomson and 
Kobori (1962) obtained compliance functions for uniformly distributed 
load over rectangular suiface area. Arnold, Brycoft and Warburton 
(1955) and Bycroft (1956) obtained the solution for rocking vibrations 
of a circular footing. 

A mixed boundary value problem, with prescnbed displacements un- 
der a rigid footing and zero tractions over the remaining surface, 
was studied by Luco and Westmann (1971), employing boundary in- 
tegral equations. Luco (1976) obtained impedance functions of a rigid 
foundation on layered viscoelastic medium using integral equations 
approach. 

Finite Element Method along with different absorbing boundaries were 
employed by several investigators to obtain solutions to dynamic soil- 
structure interaction problems. Lysmer and Kuhlemeyer (1969) de- 
veloped an absorbing boundary, called standard viscous boundary 
(SVB), consisting of discrete viscous dampers. Later they devised an- 
other boundary approximation called Rayleigh wave boundary (RWB) 
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which was recommended to be used along the vertical boundaiies 
to achieve better absorption of Rayleigh waves. They oblanu'U coni- 
phance functions for strip and axisymmetric foundations. Dasgupta 
and Kameswara Rao (1978) extended the methodology to .study plane, 
axisymmetnc and three dimensional problems. They incoipoiated 
non-linear material behaviour and non- homogeneity of the soil in the 
formulation. Lysmer and Waas (1972) developed an absorbing bound- 
ary for plane strain problems which transmits Love and Rayleigh 
waves. This boundary was applicable in frequency domain for lay- 
ered media with horizontal layers extending to infinity resting on hart! 
stratum. Kausel (1974) extended the same for axlsymmeliic problcni.s 
and obtained impedance functions for different modes. Cohen and 
Jennings (1984) developed a transmitting boundary, which they called 
paraxial boundary, by solving wave equation propagating in one di- 
rection only. This boundary was successfully employed along with 
finite elements for obtaining solutions in time domain. Underwood 
and Geers (1978) developed a boundary approximation called doubly 
asymptotic boundary element analysis, which is an exact one for very 
low and high frequencies. Higdon (1992) developed a mathematical 
approximation to boundary (Higdon boundary) absorption by employ- 
ing wave absorbers which can absorb waves impinging at any angle. 
The Higdon boundary was employed with finite difference method to 
solve wave propagation problems in layered media also. Venugopaia 
Rao and Kameswara Rao (1994) made use of infinite elemeiUs lo obtain 
the response of strip and circular footings for low frequency excita- 
tion, Zhang and Zhao (1987b), Zhang and Zhao (i987a) developed 
an infinite element approximation in frequency domain lo obtain the 
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response of strip and rigid foundations on layered media. 

The Boundary Element Method (BEM)s which can automatically in- 
clude infinite region in the modeling have been extensively used to 
study the soil-structure interaction problems. The frequency domain 
formulation was first used to obtain impedance functions of rectan- 
gular foundations resting on or embedded in a viscoelastic half-space 
by Dominguez (1978). Impedance functions of foundations resting 
on non-homogeneous and layered soils were obtained by Abascal and 
Dominguez (1986). Ahmad and Banerjee (1988) employed quadratic el- 
ements in frequency domain to obtain compliance functions Karabalis 
and Beskos (1984) and Karabalis and Beskos (1986) used a time-domain 
BEM and obtained impedance functions for rigid surface and embed- 
ded foundations The response of flexible foundations was investigated 
making use of domain methods like Finite Element (FE) and Finite 
Difference (FD) for foundation and Boundary Element (BE) for soil. 
Iguchi and Luco (1981) obtained the response of rectangular founda- 
tions on elastic half-space, KarabaUs and Beskos (1985) and Gaitanaros 
and Karabalis (1988) obtained response of surface and embedded flexi- 
ble foundations using time domain boundary element method for soil 
and FEM for the foundation, Gucunski and Peek (1993) presented 
response of circular foundations using soil stiffness matrices for lay- 
ered soils evaluated through Green's functions and stiffness matrices 
of plate using finite difference method. 

Gazetas and Tassoulas (1987b), Gazetas and Tassoulas (1987a) devel- 
oped approximate models for stiffness and damping parameters for 
horizontal and rocking vibrations based on results obtained through 
numerical computation. Wolf (1988) developed approximate models 
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by adding an additional degiee of freedom. Ho obtained stiifness, 
damping and mass parameters by curve fitting techniques to the re- 
sults obtained through exact numerical computation. Similar models 
were reported by Banos and Luco (1990) with five parameters. They 
observed that these models represent elastic half-space belter. 


1.2 Spatial and Temporal modeling 

As pointed out earlier, any realistic numerical model for dyniunic .soil 
structure interaction has to discretize the spatial as well as lime ^Ui- 
mains. Domain discretization methods (DDMs) like finite element atul 
finite difference methods etc. can handle finite size domain.s, aiul these 
methods can easily incorporate spatial variation (non-homogenily) of 
the material properties. The boundary discretization (like boundary 
element and boundary integral) methods (BDMs) can naturally adtlress 
the infinite extent of the domain. These BDMs are usually global in 
space and time. They can not easily incorporate the non-homogenity 
and non-linearity of the material. To model an infinite ilomain, by 
using a finite computational domain (in DDMs), an artificial boundary 
is to be created at a finite distance with appropriate boundary comli- 
tions to simulate smooth wave propagation ( i.e. non reflccling) at ihi' 
boundary. 

1.2.1 Absorbing boundaries 

Absorbing boundaries along with domain discretiznllon was tried .suc- 
cessfully by various researchers. The Lysmer's standard viscous and 
Rayleigh wave boundaries are local in space. The boundary scheme 
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is exact in one dimension but approximate in two and three dimen- 
sions i.e, it can fully absorb waves impinging m normal diiection only. 
The boundary can be implemented in frequency domain only, since 
the Rayleigh wave boundary description contain frequency dependent 
parameters. 

The paraxial boundary developed by Cohen and Jennings (1984) is a 
second order approximation to the waves traveling in one direction 
only. This boundary also reflects waves other than the ones impinging 
in normal direction The paraxial boundary generates unsymmetric 
damping matrices, consequently, needs special solution procedures and 
special treatment to the boundary element. This boundary was not 
found to perform better than standard viscous boundary Wolf (1988). 
An absorbing boundary, called consistent energy absorbing boundary, 
for plane strain conditions was developed by Lysmer and Waas (1972). 
They solved wave propagation of Love and Rayleigh waves in a 
horizontally layered soil resting on rigid bed rock in wave number- 
frequency domain. Solving a quadratic eigenvalue problem for wave 
numbers, corresponding to waves propagating in one direction only, for 
a given frequency, a dynamic boundary stiffness matrix was developed. 
Kausel (1974), Kausel and Roesett (1975) extended the same for an axi- 
symmetric formulation including unsymmetric modes. This boundary 
is more accurate than the previous ones, but, involves substantial 
computational effort and can be implemented in frequency domain 
only. 

The doubly asymptotic approximation boundary developed by Under- 
wood and Geers (1978) is exact only at zero and infinite frequencies. 
The boundary approximation developed by Higdon (1992) consists 
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of simultaneous appUcation of series of energy absorbing unulitious, 
le. boundary condition is expressed as a product of fust older wave 
absorbing conditions which absorb fully waves of pailicular veloeily 
impinging at specific angle. The combined effect of such boundary 
conditions is to absorb wide range of waves with dilfeienl veloiiUes. 
This boundary introduces unsymmetric damping maliicesy and con- 
sequently, numerical instabilities. These issues have to bo properly 
taken care of in the implementation using suitable procetUires like 
upwindmg. The accuracy of the scheme improves vvUh hig,her or- 
der boundary conditions and introduces more coiiiplexily aiul entails 
higher computational effort. 

Gwoli (1992) developed an absorbing boundary scheme which maps 
large domain (can be infinite also) to a fmile domain with Dirichlel to 
Neumann (DtN) map on the boundary i.e. it maps Dirichlot boundary 
condition to a Neumann boundary condition, on the arbitrary bound- 
ary created. He gave expressions for DtN map for different class of 
problems in frequency domain and time domain. In time domain, DtN 
boundary being global in time, needs to store the lime histories of the 
primary variable and its derivatives on the boundary for compulation. 

1,2.2 Boundary element methods 

The boundary element methods (exterior problems) have the capability 
to handle infinite domains. The frequency domain formulations are 
computationally competitive to the DDMs. With availability of Green's 
function for layered media (Kausel and Peek, 1982) and with combina- 
tion of FEM for Flexible foundations (GucunskI and Peek, 199.1), they 
offer a computational framework for dynamic soil-structure interaction 


1 2 Spatial and Temporal modeling 


9 


problems in the frequency domain. The time-domain counteipart (of 
BEMs) has been successfully employed to study different class of piob- 
lems like flexible foundations (Karabalis and Beskos, 1985) , vibration 
isolation (Ahmad and Al-Hussaini, 1991) and Pile foundations (Cheung, 
Tham and Lie, 1995) The BEMs suffer from following disadvantages 
They are global in space and time. Incorporation of foundation flexibil- 
ity and material nonlinearity is difficult to achieve. The time stepping 
needs nodal values to be available (both tractions and displacements) 
from initial step up to the previous step to compute the values for the 
current step ( i.e. since time stepping is global) Hence step size can not 
be changed easily and prior knowledge of proper step size is essential 
to achieve acceptable accuracy. It was established that time-domain 
formulation of BEMs did not obey causality condition and modeling 
of non-convex domains needed special treatment like subdivision into 
multiple convex domains(Von Estorff, Pais and Kausel, 1990), Due to 
the global nature of time stepping, the computational requirements 
are enormous for time-domain formulation, 

1.2.3 Temporal discretization 

Spatial discretization of the domain in dynamic soil-structure inter- 
action problems leads to a system of coupled second order ordinary 
differential equations with time as independent variable. These dif- 
ferential equations are to be integrated using numerical integration 
(time marching) schemes such as Newmark method, Runge-Kutta and 
central difference methods. 

The accuracy of the schemes depends on the step size which is gov- 
erned by the maximum eigenvalue (natural frequency) and it is to 
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be emphasized here that estimation of bounds of the eigenvalues for 
general finite element mesh is not easy. One should note lhat adap- 
tive step size methods ( i.e. the method which can predict the local 
truncation errors and suitably change the step size) are to be prefcricd 
due to their inherent higher levels of accuracy and possible savings in 
computahonal effort. 

The desirable properties of any integration scheme are unconditional 
stability, flexibility to change the step size, high order of accuracy, 
flexibility for incorporation of algorithmic damping and its spectral 
radius should be close to unity for low frequencies and liUle less than 
unity for high frequencies (Hughes, 1984). 

Newmark method is widely used in wave propagation and struc- 
tural dynamics problems due to its unconditionally stability. Dut this 
method looses its order of accuracy (second) if algorithmic damping 
is introduced. Hilber, Hughes and Taylor (1977) introduced a second 
order method which retains the order of accuracy (second) and can 
incorporate desirable levels of numerical damping. Hoff and Pahl 
(1988) derived a second order method with one free parameter which 
can incorporate numerical damping. Zeng, Wiberg, d nl, (1992) 
developed a post-processing technique to estimate local truncation er- 
rors in standard Newmark method and step size changing strategy. 
They found that this methodology Improves accuracy and reduces the 
computational effort. 
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1.3 Scope of the Present Investigation 

Based on the brief survey carried out. It can be pointed out that 
an accurate and feasible computational frame work for dynamic soil- 
slructure interaction in time-domain is not available. The domain dis- 
cretization methods along with various absorbing boundaries provide 
approximate solutions with different levels of accuracy in frequency 
domain Some of the time domain boundaries like Higdon boundary, 
provide an improvement over others but are associated with undesir- 
able properties such as generation of unsymmetric damping matrices 
The time domain boundary element and DtN methods along with 
FEM accurately represent the semi-infinite extent of the domain, but 
the computational cost is enormous due to their global nature in time. 
A computational framework which will have the desirable properties 
of the domain methods like easy incorporation of non-homogenity 
and non-linearity of the material(s), and an adaptive time integration 
scheme and with reasonable demand on computational resources, will 
be an appropriate strategy to solve dynamic soil-structure interaction 
problems in time domain. 

To achieve such strategy the present investigation uses finite element 
model for spatial discretization (for both near and far field by consid- 
ering a large mesh) circumventing the need for an energy absorbing 
boundary. Load dependent orthonormal (Lanczos) vectors (Wilson, 
Yuan and Dickens, 1979)have been utilized to transform the equation 
from original space to Lanczos space (to reduce the computational ef- 
fort). An error monitoring technique is built in the algorithm to asses 
the loss of accuracy due the transformation. An adaptive time inte- 
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gration technique has been devised by post processing the (-) method 
proposed by Hoff and Pahl (1988) with specified error tolerance. The 
response of the system is obtained in time domain by direct lime 
integration in the transformed space and requiied tcsponse quantities 
are computed in the original space. Apart from the general study of 
three dimensional foundations using 3D modelling in lime domain, 
the effect of coupled modes of vibration of has not lecievcd enough 
attenhon. The effect of flexibility in different modes of vibration need 
to be estabbshed. The effect of embeddment for different aspect ratios 
of the foundation has to be quanhfied for different modes, lifleci ol 
layering of the soil system and variation of modulus with the depth 
need greater focus. Effect of hard rock at finite deptli on Ihc responsc.s 
need to be investigated. 

The available literature on the analysis of pile foundation soil system 
when modes of vibration are coupled is meagre. The effect of pile 
cap and its embeddment has not been incorporated in the analy.sis 
which effects the responses significantly. The continuity at the pile 
soil interface was enforced approximately in the past ( i.c. continuity 
at the nodes only). 

In case of vibration isolation using open trenches the effect of depth was 
investigated. The effect of operating frequency, width and location of 
the trench need to be addressed. The mode of vibration of the source 
is also needs focus. The effect of layering and material properties of 
the fill material of the infilled trench are to be investigated. 

With the computational framework described above three dimen- 
sional block foundations have been investigated incorporating coupled 
modes, flexibility embeddment, layering etc. . Similarly the pile foun- 
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dation soil system including the effect of pile cap and its embeddment 
have been investigated for all modes of vibration. Vibration isolation 
v^ith vertical and horizontal source have been investigated The effect 
of operating frequency, width and location of the trench have been 
included in the analysis. 


1.4 Organization of the Thesis 

The work that has been carried out in this investigation is presented 
in the following manner. 

Chapter 2 describes the formulation of the present computational 
frame work - basic concepts of finite element discretization, application 
of the FEM to dynamic soil-structure interaction problems, detailed 
exposition of the concept of load dependent Lanczos vector generation, 
generation of transformation matrices, inclusion of damping, detailed 
description of the time integrahon scheme, post processing technique 
and error control strategies. A thorough validation carried out on a 
two dimensional time domain test problem has also been presented 
In Chapter 3, three dimensional embedded foundation have been 
investigated. Influence of various parameters such as aspect ratio, 
embeddment depth, flexibility of the foundation, layering, linear vari- 
ation of the modulus of the soil on the response of the foundation 
in vertical, horizontal and rotational modes are presented graphically. 
The results obtained with present investigation are compared with 
existing solutions. 

In Chapter 4 vibrations of pile foundation are analyzed. Results of 
dynamic influence coefficient for pile group effect are shown graph- 
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ically The response of the piles for transient excitnlion is coinpulc'd 
and compared with existing solutions. The effect of velocity ratio, pile 
cap etc on the response is shown graphically. A real life problem of 
a machine foundation resting on block and pile foundation has been 
analyzed. The forces transmitted by the machine consists of arbitraiy 
penodic excitahons. 

In Chapter 5 vibration isolation by open and filled trenches are inves- 
tigated. The results are compared with existing ones, The effect of 
locahon of the trench, width, depth and material pioperties of the (ill 
are investigated. Effect of layering of the soil is also investigated ami 
the amplitude reduction is presented graphically. 

Conclusions are summarized in Chapter 6. Possible extensions of the 
present work are indicated. 



Chapter 2 

MODELING DYNAMIC 
SOIL-STRUCTURE INTERACTION 


2.1 Introduction 

The dynamic soil-structure interaction deals with response of structure 
and/or ground subjected to a specified time-varying loads The loads 
could be applied directly on to the structure (arising from rotating 
machinery, wind, etc ) or through soil in the form of incident waves 
(arising out of earthquakes, explosions etc ). 

Analytical solutions for dynamic soil-structure interaction problems can 
be obtained for situations with highly idealized geometric representa- 
tion of the problem, linear elastic material behaviour, and harmonic 
loading. Since soil exist as semi-infmite half-space, numerical methods 
have to incorporate the notion of infinity in the formulation. Numeri- 
cal Methods like Integral Equation Method, Boundary Element Method 
can handle infinite domains naturally, whereas domain discretization 
methods like Finite Difference and Finite Element Methods can not 
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be appUed to semi-infinite domains perse. A way out of handling 
such infinite domains is to consider a finite (computational) domain 
for discretization with an appropriate boundary conditions (eiieigy 
absorbing) on the computational boundary, or a large domain with a 
boundary which reflects little energy. The later approach has been 
adapted in this investigation for its obvious reasons as discussed in 
the previous chapter. 

The study of Dynamic Soil-Structure interaction can be classified into 
two distinct approaches - Frequency domain approach and Time do- 
main approach. The frequency domain approach is suitable for lincai 
elashc problems with periodic excitations. The Time domain approach 
can handle all types of excitations and non linear material behaviour. 
Domain discretizahon (large size model) coupled with time integra- 
tion (time domain approach) has been utilized to analyze dynamic 
soil-structure interaction problems, in the present study. 

In this chapter the computational model utilized in this investigation 
is presented. The Modeling consist of three stages - 1. Spatial 
discretization, 2. Lanczos vector transformation and 3. Temporal dis- 
cretization. The validity of elastic continuum approach is discussed 
in Section 2 2. The Finite Element Method, which is used for spa- 
tial discretization, and its implementation is described in Section 2.3. 
The Lanczos vector transformation method, which is used reduce the 
problem size is described in Sechon 2.4. The temporal discretization 
needed to obtain the solution in time domain as implemented in this 
investigation is presented in Section 2 5. The methodology developed 
has been applied to few cases representing the range of problems to be 
addressed. The results along with published ones, wherever pos.sible, 
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nre presented in Section 2,6 for validation A brief summary of the 
conclusions came out of the this numerical modeling ate presented in 
Section 2 7 

2.2 Elastic Continuum Approach 

Soil IS a complex material whose behaviour is difficult to describe 
in simple constitutive laws. It is a particulate medium with three 
phases namely solid phase (mineral), liquid phase (usually water) 
and gaseous phase. Moreover the soil is not perfectly elastic, whose 
properties depend with depth, in-situ (initial) stress and past loading 
history, efc. Hence a unified mathematical description of its behaviour 
IS difficult The behaviour could be elastic, plastic, viscoelastic or 
thixotropic depending on loading conditions and history. However for 
small stresses, which do not induce plashc flow, an elastic description 
of the material behaviour produces satisfactory results (Barkan, 1969). 
The soil as it occurs in nature is bounded by horizontal plane at the 
top and extends to infinity in lateral directions It can be treated as 
semi-infinite body. The individual soil particles being small compared 
to the dimensions of the problem geometry, the soil is usually treated 
as a continuum. 

In the absence of any theory which takes into account the special prop- 
erties mentioned above, the elastic half-space theory is being used for 
analytical and computational modeling The elastic continuum theory 
being rigorous (not empirical) in nature, it is suitable for experimental 
verification also. The matenal properties (here soil) needed to describe 
the linear elastic behaviour — such as Young's modulus Poisson's 
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ratio V and mass density p can be determined easily. 

In this investigation soil is modeled as elastic half-space. Non ho- 
mogeneity of the material in terms of Its Young's modulus is also 
considered. 


2.3 Finite Element Method 

Many of the engineering problems encountered in practice deal with 
systems with distributed (continuously) parameters. The application of 
physics to such systems (continuua) often leads to a system of partial 
differential equations In case the equations and the domain are sim- 
ple, the solution can be obtained in a closed form (exact), often in the 
form of an infinite series, or using Fourier, Laplace transform methods 
etc. Exact solutions to general partial differential equations are difficult 
to obtain due to irregular and geometrically complicated domains and 
geometric and material nonlinearities. Consequently, various methods 
of finding suitable approximate solutions have been under continu- 
ous development. Finite Difference Method (FDM) and Variational 
methods are more popular ones. The classical variational methods 
(Ritz, Galerkin and Kantorovich methods)that are applicable to the 
whole domain, are based either on the minimization of a quadiatlc 
functional associated with the given problem or on minimization of 
error in the approximation. The difficulty in applying these methods 
lies in constructing approximating functions of the dependent variable, 
which need to satisfy the geometric boundary conditions on irregular 
domains The finite element method (FEM) which is an offshoot of the 
classical methods overcomes this difficulty and gained wide popularity 
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and acceptance (Reddy, 1984, Reddy, 1986) 

In FEM the gjven domain is subdivided into a Iinite number of small 
subdomains. The subdomains, called finite elements, are of geometrically 
simple shapes, and permit a systematic construction of the approximat- 
ing functions. Generally these approximating functions are algebraic 
polynomials developed using interpolation theory, hence called in- 
terpolation functions, and are independent of the specific boundary 
conditions and problem data Since the approximating functions aie 
defined element-wise, the accuracy of the appioximation can be im- 
proved by increasing tlie number of elements. The minimum degree 
of the polynomial used for the interpolahon functions depends on the 
order of the differenhal equation and/or the associated functional of 
the problem being solved. This in turn dictates the number of inter- 
polation points, called nodes to be identified in the element. Generally 
the nodes are placed on the boundary of the element such that they 
uniquely define the element geometry (like vertices of a triangular 
element). Additional nodes, which may be required to define the 
interpolation functions can be placed either inside or on the boundary 
of the element. The boundary nodes also facilitate in connecting the 
adjacent elements together by ensuring the primary degrees of free- 
dom (the variables that appear in essential boundary conditions) to 
be the same at the common nodes. The nodal values of the element 
uniquely defines the function values within the element and they form 
the basic unknowns of the problem. 

Though the FEM was originally formulated based on physical (intuitive) 
approach, mathematical basis was developed later . The mathematical 
approach being elegant and rigorous, is generally employed. The FEM 
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formulations are illustrated with elasticity problems. The formulation 
of finite element method can be divided into the following steps 
(Schwarz, 1988). 

Stepl The Problem domain is discretized by dividing the total do- 
main into elements. The elements in general have straight 
edges, but curved elements produce better approximation to 
the reg;ion. The region is replaced by union of approximating 
elements for computational purposes. 

Step2 For each of the elements a suitable approximation to ihe 
dependent variable which describes the problem has to be 
chosen. Rational functions of independent space variables are 
suitable. The form of approximation depends on the shape 
of the element and also on the type of problem being solved. 
These functions are evaluated using interpolation theory and 
have to satisfy continuity requirements which are generally 
obvious from purely physical considerations. The continuity 
requirements are also necessary mathematically because the 
set of approximating functions has to form an admissible cla.ss 
for either extremum principle or weighted residual method. 
If u represents displacement field of a continuous medium^ u. 
and often its derivatives (eg. for plate and beam problems) 
must be continuous across the boundary of two adjacent 
elements. Elements with approximating functions satisfying 
the conhnuity requirements are called conforming elements. 
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2.3.1 Element Models and Approximating-Functions 

The subdomains or elements makeup the domain and they approxi- 
mately represent the domain. The collection of the elements is called 
the finite element mesh. Each element in the mesh is a closed and 
non-empty one. The approximating functions as mentioned m earbei 
sections have to satisfy the continuity requirements of the variational 
statement of the problem. The approximating functions/ depending 
upon their degree, could be linear, quadratic, cubic polynomial func- 
tions of the dependent variable. If u' represents an approximation to 
the displacement field defined in the domain of the element, it can be 
expressed by a degiee polynomial of the independent variable 

T as 


= tti -F a2X -f a^x'^ -t- • + apX^ ^ (2.1) 

The coefficients aj, 02 etc can be evaluated using the function values at 
the nodes ( t.e, nodal variables) of the element. Then the approximating 
function u® can be expressed as a linear combination of basts functions 
with nodal variables as coefficients as follows. 

%• =.f^u‘N;(,x) ( 22 ) 

J=1 

where denotes the nodal values of the dependent variable u®, and 
Ni are algebraic functions of x and are of degree p - 1 . The above 
equation must be valid for any nodal variable of node i; the basis 
functions N^(a;) must possess the interpolation characteristic that at 
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Convergence requirements : 

In any numerical formulation, solutions obtained must converge or tend to the exact 
solutions of the problem. In displacement formulation the approximating functions 
liave to satisfy certain requirements to ensure the convergence (Desai and Abel, 1972). 

1 . The displacement models must include the rigid body modes. This can be 
achieved by liaving constant tenns in the model ( like at in equation 2.5). 

2. The displacement models should include the constant sUain states of the element. 
As element size approaches zero, the element strain approaches a constant value. 
To incorporate tins behaviour the model must iiave terms associated with constant 

i 

strain ( like ai in equation 2 5). 

3. The displacement models must be continuous within the elements. The continuity 
requirement can be met by choosbg polynomial models for approximating 
functions. 

Elements satisfying both second and tliird conditions are called complete elements. » 
Elements satbfying patch test will converge (Reddy 1986). 
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Element Models 

The elements can be classified based on the degree of llie approxinuU- 
ing (interpolation) functions. They could be linear, quadralic, cubic dc. 
The interpolation functions can be consliucled using Lagianj>e inter- 
polation, Hermite interpolation, etc. The piinciples would be illustrated 
with a one-dimensional elements and the procedures for constructing 
the interpolation functions for two and three dimensional cleinenls 
will be indicated. 

1-dimensional Elements 

A typical one-dimensional element family is shown in higure 2.1. I, el 
■u® be dependent variable and x the independent variable, hor linear 

t 

element of length h, can be approximated as 



Qimdrnllc LMctnciils 

(a) Global ooorUmalc sy,le„, 

Fiffure 2 1’ Typical Oim -diiTieiisinnal [’j('nip|ils 
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((' «! -I ^2.5) 

Constants O) and can be evaluated utilizing the nodal variables of 
</' al nodes 1 and 2 in the above equation, te 

:v -- 0, (f ~ v\ = d] 

-- h, II ~~V2-- tii -b02h (2 6 ) 

Solving (lie above equation for coefficients ai and the ( 2 . 5 ) can be 
wiillcn ns 


V*' N)Ui -)- N-jfU-j 

Putting the same in matrix notation 


( 2 . 7 ) 




«i 

^2 




wlierc 


( 2 . 8 ) 


I- x/h and N2 = x/h ( 2 . 9 ) 

Ni and N2 are called shape functions because they define the shape 
of the variation of 11®. L^®J is a row vector of shape funedons 
containing N\f N2 and {11®} is a vector of nodal variables ui, 112. As 
already mentioned in the earlier section, the shape functions can be 
constructed making use of interpolation theory (such as Lagrange s). 
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For example for a three noded element, the appioximalins function 
can be written as 


= L^i, A^3,A^2j \ 


U\ 


(2.10) 


Now the shape functions N^, N-i, Nz are quadratic functions of the 
independent variable x, and can be constructed using I.agranj’o's 
interpolation as follows 


'v.=n 




(a:, — Xj) 


( 2 . 11 ) 


Where J] indicates product of all possible terms and n denotes number 
of nodes of the element. 


2 and 3-dimensional Elements 

Shape functions of two and three dimensional elements can be con- 
structed as products of corresponding shape funclion.s of one the 
dimensional elements in each direction. 
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Linear LlcnKni^ 



QuaJraili ICIcinciKs 


(a) Global coordinate sy&icm 


(b) Nalural coordmalc systen'i 


Figure 2 2 Typical Two-dimensional Elements 

Typical rectangular elements of linear and quadratic interpolation are 
shown in Figure 2.2. The shape functions Ni to Ni for the linear 
element can be constructed from its 1-dimensional counterpart as 
follows: 

Ni(x,y) = Ni{x)*Ni{y) 

= ^2ix) * Ni{y) 

-^3 (a;, y) = N2(x) * N^iy) 

t/) = Ni{x) ♦ N^iy) (2.12) 


Quadratic element shape functions can be constructed in a similar 
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fashion. Tho serendipity elements/ which are special laniily ol clcniL'iils 
(Figure 2.2) will have nodes only on the boundary of the elemonl. 
Construction of shape functions for the higher order (2 and above ) 
serendipity elements need special treatment (sec HugheS/ 1987 for 
details). 

Natural coordinates and isoparametric elements 

To facilitate generalization of computation of the shape funelions and 
element matrices a local (element) coordinate system will be bene- 
ficial. A natural coordinate system is a local coordinate .system in 
which a point is specified in terms of dimensionless numbers whose 
magnitude does not exceed unity. A natural coordinate system for 
one-dimensional element is shown in Figure 2.1. The relationship 
between natural coordinate ( and the physical coordinate .r can be 
written as 


2x - {X 2 + Xi) 

X2 - Xl 


(2.13) 


As earlier pointed out the accuracy depends on the approximation of 
the domain by individual elements. To represent curved bouiularies a 
coordinate mapping in which an element described in natural coordi- 
nate system is mapped on to a curved element described in physical 
(global) coordinate system. If the interpolation functions used for de- 
pendent variable approximation and the functions used for coordinate 
mapping are same then the elements are called isoparametiic eiemenls. 
The coordinate mapping in one dimension can be expressed as 
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•T = [NiNil 


•'ll 


'12 


[WH-.} 


(2.14) 


Now the interpolation (mapping) functions A^i and can be rewritten 
in terms of natural coordinate ^ as 


= i(: - 0 and Wj = i(l + 0 (2 15) 

Similar expressions can be developed for 2 and 3-dimensional elements 
analogously as follows. Typical isoparametric elements m 2-dimensions 
are shown in Figure 2.2. The coordinate mapping m two dimension 
can be expressed as 


a;i 

Vi 



Ni 0 

N2 

0 

Ns 

0 

N4 

0 

\yf 

0 Ni 

0 

N2 

0 

Ns 

0 

N4 

J 


2/2 

■' 1‘3 


2/3 


(2,16) 


X4, 


/ 


The above equation can be generalized for 7i-dimensional space as 


X = [N,] { X, } (2.17) 

where x and Xj are the vectors of independent variables, and coordi- 
nates for the node x respectively; [N^] is a diagonal sub matrix of order 
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equal to the number of space dimensions. Kach of ita entry equals to 
mapping function of the node Expressing the s.iinc ni.itheinatieally 


■x. — \xyz 



[N.l = 


N, 0 0 

0 0 

0 0 Af, 


0 

0 

0 


[ 0 0 0 '•< iV, 

rp 

X, = 1^‘t 2 /, . . ] 


(2.1H) 


The interpolation (mapping) functions N\ to A'^i of the bilinear (lineai 
in each direction) quadrilateral element (I'igurc 2.2) can be written 
from it's one dimensional counter part (linear element) in natural 
coordinates as follows: 


^3 = ^(1 + C)(l + 7/) 

= -<)(' + '() ( 2 . 1 ')) 

Since approximating functions are same as the mapping functions In 
isoparametric element scheme, the dependent variables can bo written 
in terms of their nodal values as 
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here u and ii, are arrays of dependent variables and it's nodal values 
lespectively 

2.3.2 Formulation 

The basic formulation can be arrived at using different methods Some 
of the popular methods are — 1. Ritz foimulation and 2. Weighted 
Residual formulations. In Ritz formulation a suitable functional is 
constructed from the physics of the problem, wheie as the differen- 
tial equation governing the behaviour is used in weighted residual 
formulation. 

For elasticity problems a suitable functional can be obtained using 
principles such as Principle of Minimum Potential energy (displacement 
formulation), Principle of Minimum Complimentary Energy (Stress 
formulation), Hamilton's Principle (dynamic problems) and Hellinger - 
Resissner Principle (Mixed formulation). The functionals so constructed 
using the above principles can be used with Ritz type foimulation 
Alternatively Differenhal equations governing the behaviour can be 
used with weighted Residual Formulation. 

Ritz Formulation 

In Ritz Formulation a functional with stationary properties is con- 
structed by applying ohe of the above mentioned pnnciples. The 
methodology is illustrated applying the principles for a static elastic- 
ity problem in the following paragraphs. Let 11 denote the potential 
energy of an elastic body, it can be expressed as 
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(2 21 ) 


The potential energy of the whole body can be expressed as the sum 
of the potential energies of the individual elements nrl 


n = 





( 2 . 22 ) 


where u, e, a are vectors of displacements/ strains and stiesses; / and 
T are the vectors of body forces and surface tractions lespeclively. il 
denotes the domain of the body and T is the boundary of the body. 
The variables with superscript e refers to the corresponding elemental 
quantities. 

Using the linear constitutive relations/ stresses can be expre.sseil as 


a = (£;]€ anda^ = [E]6‘ (2.23) 

where [E] denotes elasticity matrix. 

Assuming small strain theory, strains can be defined in teim.s of dis- 
placements using differential operator (T) as 


^ (2.24) 

The differential operator [L] depends on the nature of the problem. 
For 2-Dimensional problems, it can be written as 
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0 




0 


Ot 


A jO 

I Oy 8r \ 


(2,25) 


Similar relations hold for elemental quantities also, making use of the 
approximation (Equation 2 8) of the element displacements; element 
strains can be expressed as 


= [Ly = [L] [fVJ = [5] y} (2 26) 

where matrix [B] relates element strains to nodal displacement vari- 
ables. 

Substituting the expressions for stresses and strains in the (2.22) results 


nel 

n=E 


1 / 

Jn‘ Jv‘ 


(2.27) 


where 0.^ and P denote the domain and boundary of the element. 
Taking variation of both sides and setting it to zero (stationary condi- 
tion) we get the algebraic equation governing the behaviour in terms 
of the nodal variables. 


nel 


6n = 


l=l 


- / ~ I 

v/m •/r* 


= 0 (2,28) 
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Since the variation 


(5u'} is arbitrary the expression inside llie paienlhe- 


sis should vanish, resulting in an algebraic equation m nodal variables, 
i e. 


= (P) 

where 


(2.29) 


nel p 

IK] = ^ I/r] and [/r] = J l/3r'[K)(/i)fh) (2.20) 

{-1 

and 

net p 

{/>} = r {/) and {p”) = / - / \N\' {■,),!« (Ml) 

" Jn« -h'" 

[/<■'] and {p®} are called elevient stiffness matrix and clement load vecloi 
respectively; [K\ and (P) ate called sU((v\ess and load vector 
respectively. The summation sign in equations (2,30, 2.31) indicates 
assembly of the corresponding element quantities. 

Weighted residual formulation 

Weighted residual formulation can be used if construction of a func- 
tional is not obvious or difficult as in the case of nonliiuMi' probli'in.s. 
This formulation makes use of the differential equations governing the 
behaviour. The approximating (unctions of the dependent variable 
(solution) are constructed as explained in the earlier section. This 
approximate solution is substituted in the differential ecjualion and 
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the residue is minimized over the domain (each element) Depending 
upon the choice of weighting functions foi minimization of the lesidue, 
different schemes are developed. In Galerkin's method interpolation 
functions are used as weighhng functions This method generates 
symmetric stiffness matrices for self adjoint diffeienlial opeiatois Let 
the equation governing the behaviour of the problem in the domain 
with the boundary conditions on the boundary F, be 


C{u) = 0 in domain Cl (2 32) 

C('u) = 0 on the boundary F (2 33) 

Where L and C are differential operators. 

The residue R over the whole domain fi is the sum of the residues 
of each element of domain Cl^ which is obtained by substituting the 
approximate solution defined over each of the elements in the above 
equation, t.e 


nel ncl 

R = Y, R^=:YlC{u^) (2.34) 

1=1 t=i 

Approximate solution can be obtained by minimizing the residue with 
respect to a set of weighting function. The weighting functions are 
chosen to be shape functions in Galerkin's method. Now the unknown 
nodal variables appearing m the approximate solutions can be obtained 
by solving the resulting algebraic equations, 



'^Rdn = 0 


(2.35) 
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This can be expressed in terms of element qiuin lilies as 





(2.J0) 


The continuity requirements of the shape function (approxiinatinp, 
functions) depend on the order of the operator C. This can be 
reduced by carrying out the integration by paits. The above et| nation 
produces a set of simultaneous equations in unknown nodal variables. 


2.3.3 Application of FEM to Dynamic problems 

In this section the implementation of the Finite lilemenl Method for 
spatial discretization is discussed. A displacement type lornuilation 
uses nodal displacements as primary variables. Varialioiuil 

formulation has been utilized to develop the procedure foi spatial 
discretization. Hamilton's principle has been used to construct the 
functional. The principle is stated as follows (Desai and Abel, 1V72) 
Among all possible time histones of displacement conjigitmlions which satisfy 
compatibility and the constraints or kinematic bounilnri/ conditions and lohk/t 
also satisfy conditions at time ti and k, the history which is ihc nclual 
solution makes the lagrangian functional a minimum. 

This can be stated as 


6 



Cdt = 0 


(2.37) 


where 6 indicates the variation of the integral and C is the l.agiangian 
which can be expressed as 


2 3 Finite Element Method 


37 


C = T -P-W (2.38) 

where T, P, and W aie kinetic energy, strain energy and potential 
energy of the applied loads Accordingly 

£ = i pv — o\ dv -V f‘v? fdv-\- f rds (2 39) 

2 Jo. Jo Jy 

where p is mass density of the material, dots denote the time derivative 
of the variables. 

Lagrangian L of the whole domain can be written as sum of the 
contributions over each of the elements, i e. 



(2.40) 


Substituting for the approximate displacement vector defined over 
each of the element and stress strain relations in the above equation 
results 




(2 41) 


Hereafter the summation symbol is omitted for the sake of conciseness 
and wherever elemental quantities appear, usual summation (assembly) 
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over the elements is implied. Substituting the above expre.vslon foi 
Lagrangian in (2.37) results 


- I IWJ’’{/'}dv - / 1,iVJ’'{t‘).I.s) ,11 II (2.-12) 
yu' yr* 


Integrating the first term 



Jn^ J 

(it = 


rt2 

Jl\ 



tit (2.43) 


substituting this expression in the (2.42) results 


h 


{6u^} f 

/^_(pW’'[tvj(ir) + |n|'-ifci|N|(a"}^</ii 

- / ltVJ'''{/')tlv- / LWJ'''(r"}-/.v 
Jq.‘ J\'< 


til 0 (3.44) 


The first term in the right hand side (r.h.s) vanishes because becaiise the 
displacement configurations should satisfy the boundary condition.s at 
times ti, and t2 (Hamilton's principle) ie. the vaiintion (Z)!!''} should 
be equal to zero at times fi, and t2. Rest of the expres.sioii inside the 
parenthesis should also vanish because the variation i.s arbitrary. The 
resulting equation can be written as 
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where is element mass matrix and is given by 


(2.45) 


/ p[AfJ^'[AJdv (2,46) 

Jnc 

Assembling the element matrices results in a system of differential 
equations (equations of motion) governing the behavioui of the body 
as follows 


{M]{U} + [iq{V} = {P} (2.47) 

Here the vectors t/, and P are time dependent. 

Mass Matrices 

The mass matrix defined by the above equation is called consistent mass 
matrix because the same shape functions [A/J are used to generate the 
stiffness matrix Often a diagonal form of the mass matrix is adopted 
assuming that the mass of the body is concentrated (lumped) only at 
the nodes of the body i.e. the continuum is discretized a§ particle 
lumps. This form of the matrix is diagonal and is called lumped mass 
matrix. The procedure is called adhoc lumping. Numerical integration 
(special) rules which produce a diagonal consistent mass matrix are 
possible for simple elements only. Usually it is recommended that con- 
sistent mass matrices be used for natural frequency computations and 
lumped mass matrices be used for wave propagation problems (Cook, 
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1974) Diagonal mass matrices are compul.Uionally aclv.uilageous .mil 
require less storage. In this investigation a systmnhc lumpms^ mol hod 
IS adopted which performs better than both adhoc lumping and con- 
sistent mass matrices In this method only diagonal onliics of the 
consistent mass matrix are generated and they aic scaled so that mass 
of the element is conserved. 

Damping 

Damping in dynamic soil-sliucture interaction is present in two form.s 
- 1. Radiation damping and 2. Material damping. The radialion 
damping is due to the infinite extent of soil in which energy is lost in 
the form of outgoing waves. Analytical formulations can incorporate 
the infinite extent of the soil, but computational methods which work 
with finite domains like FEM can not incorporate the same. As ahoady 
mentioned the concept of infinite extent can not be incorporated in FliM 
formulation and hence A Large Mesh is considered lor discreli/alion. 
The energy at large distances is negligibly small and the reflections 
do not influence the accuracy of the solution.s. The material daiu]nng 
is due to hysteresis in the material. The Irentmcnt of the material 
damping in computational analyzes Is often specified in tt'rms of 
equivalent viscous damping. This is introduced by rnean.s of specified 
fractions of critical damping, or Rayleigh or proportional damping in 
which damping matrix is specified as a linear combination of ma.ss 
and stiffness matrices. The former description is advantageous in 
frequency domain analysis, whereas the later form Ls used In lime 
domain analysis. The Rayleigh damping is utilized in the pre.sent 
investigation. The damping matrix (C'] can be expressed as 
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\C] = cx\K]'\- P\M\ (2 48) 

where a, and /? are called the stiffness and mass pioportional damping 
constants respectively. Incorporating the damping defined as above 
the equation of motion becomes 


\M]{U} + (C]{L/} + [K]{U} = {P} (2.49) 

The above equation which describes the behaviour of the dynamical 
system is transformed to Lanczos space as explained in the next section. 
The transformed equations are used in time stepping scheme. The 
time histories of the generalized coordinates (of transformed domain) 
are projected on to the original space to obtain the displacement 
histories etc. Hereafter the parenthesis around a variable [] indicating 
they represent a matrix will be omitted for conciseness. 

2.4 Lanczos Method 

2.4.1 Introduction 

Lanczos method is originally devised to extract some (not all) eigen 
values and eigen vectors of an eigen system It involves generation of 
orthonormal sequence of vectors known as Lanczos vectors iteratively 
from Krylov subspace. At each step the current vector is normalized 
with two previous vectors and build a tridiagonal matrix whose eigen 
values and eigen vectors are Rayleigh-Ritz approximations of the 
original system. This method (theoretically) terminates after n (order 
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of the system) iterations. But, in practice, due to Unite pieeision of 
numerical computation, the Lanezos vectors loose their oilhop,on.ility 
and may need re— orthogonalization frecjuenlly (see laij-^e, L97(i). I he 
same procedure is extended to generalized eigen value piohleins. 
The loss of orthogonality of Lanezos vectors is investigated by Laige 
(1976). Different schemes -- total re-orthogonalizalion (WiKson, Yuan 
and Dickens, 1979), selective re~orthogonalizalion (Parlelt and Scott, 
1979), and partial re-orthogonalization (Simon, 19H4), are suggesleil to 
overcome the problem. 

The Lanezos vectors which are orthonormal, can also be used to 
construct a coordinate transformation matrix to be used in Ihe .solu- 
tion of dynamic equilibrium equations. Recently Wilson, Yuan and 
Dickens (1979) and Ibrahimbegovic, Harn, et uJ. (1990) developetl a 
methodology to generate a set of orthogonal (Lanczo.s) vectors (or use 
in coordinate transformation as an alternalive to conventional mode 
super position method. Nour-Omid and Clough (1984) sluivved that 
Lanezos vectors can be used to solve slructuial dynamic problems in 
time-domain. The basic methodology as implemented in this invesU- 
gation in generating Lanezos vectors is dc.scribc(.l heie brietly. 


2.4.2 Lanezos algorithm 

The Lanezos vectors are generated spanning the Krylov sub space de- 
fined by the sequence of vectors r, {[<' LV/)‘ r , (A' L\/)'* i\ {K ' Aff r 
, . (A'-LV/)‘”r, given the pair of matrices A and Af and r being 

an arbitrary vector. The sequence of vectors convcrgo.s to an eigen 
vector corresponding to the smallest eigen value of the generali/ed 
eigen value problem defined by the equation. 
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(A ^ XM)z = 0 (2.50) 

111 dynamic response analysis the oithonormal sequence of vectors 
from Krylov subspace will be utilized for coordinate tiansformation. 
Different schemes (see Cullum and Willoughby, 1985) are available foi 
generation of Lanczos vectors for eigen vector computation. In Lanczos 
algoiithm, Gram-Schmidt orthogonalization (Cullum and Willoughby, 
1985)will be employed at each step to normalize the current vector 
with respect to matrix M and the two previous vectors of the Krylov 
sequence The result will be a set of m orthonormal vectors 
Let (gi,< 72 , ,(}}) he j Lanczos vectors generated and vector is to 

be found Then can be calculated by hrst computing a preliminary 
vector from the BCrylov sequence 


fj = K~^Mqj (2.51) 

Now the preliminary vector will be normalized with two preceding 
vectors qj- 2 i gj~i as 


+ 03<l3-i + + ■ • (2.52) 

where is pure vector orthogonal to two previous vectors and 
Oj, f3j and 7 , are the amplitudes of the previous vector contained 
in fy The coefficients can be evaluated using the orthonormality of 
Lanczos vectors. Premultiplying the above equation by gjM , and 
noting the M orthonormality of the vectors, except first term all other 
terms on the right hand side vanish, resulting in 
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aj = qjMrj 

Similarly /?, can be evaluated as 


(2.ri3) 


P] ] 

Replacing fj by it expression from (2.52) we get 


(2.54) 


p, = = rf. yM(jj (2 55) 

Expanding 7"j_i in terms of its pure vector I'j i sub.sliluling in llie 
transpose of (2.55) pj becomes 


Pj = Mr + aj-iqjMqj-i + pj..iqjM<ij ^ I- • • • (2.56) 

It is obvious that that all terms except the first vanish on the right 
hand side of the above equation. Now i/j is the vector ohlaincd by 
normalizing rj_i with pj, ie. 

1 

% = ^0-1 (2.57) 

Substituting this expression in (2.56) results 
1 

P] = (2.5H) 

Then the value of will be 


P] = r,_iiV/r_,_i 


(2.59) 
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In a similar procedure it can be shown that - 7 ^ and rest of the coefficient 
of the terms in (2 52) can be shown to zero 

In summary, the process of generating the Lanczos vectors can be 
expressed as follows: 

Choose 9 o = 0/ a starting vector tq and fit = . then, 

For j ~ 1,2, ,m compute Qj 


Tj-l 

(2 60) 

II 

1 

(2 61) 

Tj — ^j(}] 

(2 62) 

aj = f jMgj 

(2 63) 


(2.64) 


2.4.3 Starting vector 

The starting vector rp, in general can be chosen arbitrarily. If the 
starting vector is related to loading amplitudes we may expect that the 
modes which are contributing to the response will be included. Thus 
the starting vector is taken to be static displacement vector given by 


ro = K-^P (2.65) 

where P is the vector of loading amplitudes. Inclusion of the static 
displacement vector (Wilson, Yuan and Dickens, 1979) avoids any possi- 
ble need of applying static correction, because the static displacements 
are included in the coordinate transformation itself. More over for 
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small number of Lanczos steps, the basis vectors will be I^ayleij\h Kil/, 
approximations of scattered eigen vectors non orthogonal to loading 
amplitudes in the Krylov sequence. The coordinate tran.sfoiinalion 
matrix will contain approximations of eigen vectors ol hij’her modes 
excited by the loading, resulting in a small number of basis vectors. 


2.4.4 Loss of orthogonality of Laaczos vectors 

As mentioned in earlier section the Lanczos vectors loose oilhogonality 
due to the finite precision aiithmetic. The simple I’aige style alp.oiithm 
(Paige, 1976) without re-orthogonalizalion may need large luimber 
of vectors (m » n) to isolate eigen vectors. Ilowevcu' the use of 
Lanczos vectors can only be pursued in coordinate Iransfornuition if 
and only if few orthonormal vectors could produce results aceeplable 
of accuracy. Hence re-orthogonalizalion of Lanc/os vectors is es.senlial 
in the application of coordinate transformation process. 

The orthogonality of the Lanczos vectors can be maintained by re - 
orthogonalization. These methods can be categoiized as - 1. (’om- 
plete re-orthogonalization 2. Limited re-orthogonalization. In com- 
plete re-orthogonalization, the current vector Is normalized with all 
previous vectors (Wilson, Yuan and Dicl<en.s, 1979). Tlil.s proce.s.s iieed.s 
large core memory or secondary storage to store all previous vectors. 
Consequently the procedure is both I/O intensive and needs large 
computational resources. In the latter choice the re-orthogonali/iUion 
IS applied when found necessary. In limited re-orthogonalizalion the 
loss of orthogonality is monitored continuously (at the etui of each 
step). Two variants of limited rc-orthogonalization namely selective 
and partial re-orthogonalizatlon schemes, which arc useful in coordi- 
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nate transfoimation process are presented. 

Selective re-orthogonalization 

This method is suggested by Parlett and Scott (1979), wherein, re- 
orthogonalization is done against a few selected (not all previous) 
vectors when found necessary It is suggested by Paige (1980) that loss 
of orthogonality is due to presence of components of converged Ritz 
vectors of the system in the Lanczos vectors. This fact is made use 
of in selective re-orthogonahzation. The current vector is normalized 
with respect to all previously computed converged Ritz vectors. These 
Ritz vectors will be small in number, and they may be stored in core 
memory. A refined methodology is suggested Parlett and Nour-Omid 
(1985) to monitor the loss of orthogonality. The computation of Ritz 
vectors is carried out periodically rather than at the end of each step. 

Partial re-orthogonalization 

The aim of all re-orthogonalization schemes is to prevent loss of 
orthogonality i.e to maintain a certain level of orthogonality among 
Lanczos vectors. We define level of orthogonality k among the the 
Lanczos vectors at step as 

K— max (2.66) 

i<r,<j-i ^ 

The full re-orthogonalization of Qj against all previous vectors maintains 
level of orthogonality k a round of level e. However numerical results 
suggest (Parlett and Scott, 1979, Simon, 1984) that semi- orthogonality 
{.e. K = \/l, among the Lanczos vectors is sufficient to prevent spurious 
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eigen values. This scheme suggested by Simon (19Hd), uminlain.s .seiui- 
orthogonaUty among the Lanczos veclois. The monilorin|- ol In.s.s 
of orthogonality is done by updating an army consisting of levels 
of orthogonality. These levels of orthogonality are not eompulecl 
explicitly rather they are computed using a simple recurrence relation 
suggested by Simon (1984). Let u 's refer to the level of oilhogonalily 
between the Lanczos vectors referred by its subsciipls, Then 


tOf^ n =1 for k — 1,.. ,j 
tefc.il - 1 ~qUk~i for 

i3j+lV'jl-l,fc=i9fcHt*^j,fc+lT{0!fc~«j)Wj,fc-l /lfcW,,c 1 {IjlO, 1,^ I 

(2.67) 

where is a random number to account for roundolf eriors. U 
is suggested by Simon (1984) that V-’;, fc ^'>>1 be lakun to be I’liiKlom 
number with 0- mean and c - standard deviation, lluie it .should be 
noted that the oj 's are estimated or computed levels of ovlliogonalily 
but found to be very close to the actual ones (Simon, 1984), 

At stage, update the recurrence relation for W/m ami fnul the 
vectors whose levels of orthogonality exceeded the si’iul-orlhogoniilily 
level. Le. 


for A: = 1, 2, . j find /s such that |[wj , i. nil > >/< (2.68) 

The vectors which satisfy the above ccjuation arc rO' oilhogonali/.ud 
with current and next vector and corresponding w.v are ro.sid to f. 

It IS found that both Selective and Partial reorthogonali/allons rtM|uiu' 
comparable computational effort (Simon, 1984). I'he Selective re- 
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orthogonalization is more suitable for eigenvalue computations since 
Ritz vectors are computed periodically for re-oithogonalization Wheie 
as Partial re-oitliogonalization schemes is suitable in case Lanczos 
vectors are intended to be used for solution of simultaneous equations 
and coordinate transformahon Here, m this investigation Partial re- 
orthogonalization scheme is chosen because of its simplicity 

2.4.5 Number of Vectors 

The number of vectors required to obtain the results with desired 
accuracy can be obtained as follows According Bayo and Wilson (1984) 
and Coutinho, Landau, ef ah (1987) a measure of error introduced 
in the the loading function in the transformed space can be used to 
gauge the accuracy of representation. It is expected that lower the 
error higher is the accuracy. Let the transformation matrix be be 
defined by the relation 


U = <pX (2.69) 

where X is the generalized displacement vector in the new coordinates 
space and 0 = [<^i, <^ 2 ) • m < nis the coordinate transformation 

matrix. The dynamic response C/(t) is now approximated by a linear 
combination of generalized displacements X(f). Substituting for U in 
(2.49), we get 

MX + CX + KX = P{t) (2.70) 

where M, C, K and are generalized mass, damping, and stiffness 
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matrices and Pit) is the generalized load voLtor lespci lively .uu! mo 
given by the following equations. 


II 

-8- 

(2 71) 

C — (f}^C(j> 

(2.72) 

K = 4>^'K^ 

(2.7:i) 

(1 

-8. 

(2.7-1) 


Premultiplying the (2.71) by M(f) and consiclciing Iho .spatial vaiialuni 
only, the load representation P,n in original space inoorpoiali'd in Iho 
analysis can be expiessed as 


P,„ = ^M<p,p, (2.75) 

J = 1 

where pj is jtk participation factor and is givL*n by 

= (2.70) 

Now the error E in loading representation can bo oxpros.socl as 


E^P-P„, 


(2.77) 


It is convenient to express the error in suitable norm rulalivo In Iho 
original load as 


e = !ra! 

IIP'^PII 


(2.7H) 
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The relative error norm e will lie between 1 (for ??j = 0 re 0 vectors 
are used) and 0 (for ??? = 7i le. all vectors are used) A suitable value of 
error tolerance (relative viz in terms of percentage) may be specified 
and the Lanczos vectors are geneiated till the relative error norm e 
falls below the error tolerance of the load vectors. 

2.4.6 Coordinate transformation matrix 

Coordinate transformation matrix can be deiived as follows: 

• Triangularize the stiffness matrix K, i.e 

K = L^DL (2.79) 

• choose static displacement vector as ro r as explained in the 
previous section 

• generate m Lanczos vectors which are orthonormal (apply partial 
re-orthogonalization where necessary). Check whether e relative 
error norm of load representation is less than error tolerance 

• Project the stiffness matrix in to Lanczos subspace as 

K„,:=Q^'KQ, g = (91,92,. >Qm] . (2 80) 

• Solve the eigen problem in Lanczos subspace 

KmV = y = bn J'z. . ■,yn.] ( 2 . 81 ) 

• Compute coordinate transformation matrix 

= gy (2.82) 

asi>J’rr?AL uks.viiUV 

I. 1. T., iCADPUU 

125692 
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2.5 Direct Integration 

2.5.1 Introduction 

Dynamic analysis of engineering systems reqiiiio temporal aiul otten 
spatial discretization (for continuous systems). The spatial disci eti/a- 
tion results in a system of coupled ditfeiential etiualions such as (2.^19), 
describing the physical behaviour. These eipiations, in jirinciple, can 
be solved by direct integration. The direct integration method will be 
efficient m case the excitation is either transient or peiiodic consisting, 
of combination of several frequencies. I’or nonlinear ilynanik* prohUuns 
the direct integration methods are the only tools readily a}'plicable. Di- 
rect integration schemes can be broadly divided inlo Iwo ealegguies 
1. Single step methods and 2. Multi-step melhotls. Sin)',le step nu'th- 
ods have a number of advantages over mulli-slep methods. They are 
self starting and amenable to time step modification. Tlu’se methods 
are more convenient for nonlinear systems. 

Linear multi-step (including single step) methods can be classified 
into two broad categories with reference to stability. As pei ITililqiiisl 
(1963) 

• An explicit algorithm which is unconditionally slaldi' does not 
exist. 

• Unconditionally stable algorithm of order greater Itum two does 
not exist. 


Thus one category of algorithms are implieil with uncoiulitional sta- 
bility of order one and two. Newmark methods (order two) fall uruim- 
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this category Zienkiewicz, Wood, ct al (1984) have developed a uni- 
fied set of single step integration methods which are unconditionally 
stable These algorithms are derived by weighted residual appioach 
In structural dynamic analysis computations damping out higher fre- 
quencies(called numerical damping)is a desirable property m addition 
to unconditional stability The spectral radius of higher frequencies 
should be less than 1 0 , whereas for low frequencies it should be 
close to 1 0 . This can be achieved by making 7 > 0 5 in Newmark 
methods, but order of accuracy falls below two, Hilber, Hughes and 
Taylor (1977) developed second order accurate method which they 
call a. method with an additional parameter a. This method has the 
property of controllable damping (specified in terms of a). Another 
method (called © method) based on weighted residual technique is 
presented by Hoff and Pahl (1988) which has better order of accuracy 
and controllable numerical damping. The © method has only one free 
parameter 6, which controls numerical damping. 

The other category falls under conditionally stable multi-step algo- 
rithms of order greater than two. Because of conditional stability and 
increase of the order of matrices by a factor of two, this class of algo- 
rithms are not investigated for use in structural dynamics. The higher 
order algorithms, because of higher accuracy, are viable alternatives if 
stability is ensured, produce more accurate solutions with less compu- 
tational effort. However, as mentioned earlier, the desirable property 
of numerical damping, can not be incorporated in higher order meth- 
ods hke 4^^* order Runge-Kutta method and Adam-moulton methods 
(Stoer and Bulrsch, 1980). 
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2.5.2 0 metliod: 

The 0 method of Hoff and Pahl (1988) which ha.s been ii.seii m Ihi.s 
investigation, is desciibed biiefly. 

The dynamic equilibrium may be described by the (I'.qua- 

tion 2.49) or by the transformed (2.70) with the given inilial conditions. 
If Xn, v„ , a„ and P„ are the displacement, velocity, acieler.Uion and 
force vectors at nth step respectively; then the corresponding quantities 
at (n + 1)^'^ step (with step size li) can be obtained as 


^R+l 

Vn+l 


^»+l 

Mjn^a 


(1,1 H‘ 6(1 

Vn + /ia„ + (15- 0)h6a 
x'n + hVn + 0.5h\ + 

dOf 

Pn f 1 M(l)x C/Hii (1 /v .1',, I 1 


where 

iLr„. = M + (15-0i)/iC + -~-/i'‘*A- 
Pni-l = Pn + t^o(P,i (-1 “ P„) 

^nhl = Un+0ik(t„ 

•Tn+i = x„ + 0ihv„-\-{).6h\x 


(2.Hd) 
(2 84) 
(2.Hh) 
( 2 . 86 ) 


with 0 95 < 01 < 1.0 

The matrices K, C, M are the stiffness, damping and mass iniUiiees of 

the system. P„ is the load vector at idli step. The.se cpianlltles refer 
to the transformed system. 

Since this is a second order method, the step size .tdjustnicnl cm bn 
incorporated using Richardson extrapolation (Stoer and Dulrsch, fJtil)). 
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That is, compute the solution at time i + h with step size and A/2 
Then the difference of the solutions gives an estimate of error as 


= (287) 

Where eir is a measure of local truncation error, a:/, and i /,/2 are 
displacement vectors at time t+h with step sizes h and h/2 respectively 
The symbol f| (j here denotes a maximum norm. 

2.5.3 0 Metliod with post-processing 

This method differs from the earlier one only in computation of local 
error. The post-processing of solution is done at the end of each 
step to estimate the error and consequently the step size required for 
prescribed error tolerance In this method acceleration is assumed to 
be constant in each step If a continuous ( viz. linear) valued function 
of time is utilized for approximation of acceleration, more accurate 
results could be obtained (see Zeng, Wiberg, et ah, 1992) If a* is the 
acceleration at any time t during the nth step, it can be expressed as 


a* = an + ~5a ( 2 . 88 ) 

where tn <t< fn+i and li = t„+i - tn 

Now the error in acceleration during the step can be written as 


a* — a 


(2.89) 
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where a is the average acceleration vector dining the .iiul can 

be expressed as 


a= 1- 


2dl 




( 2 /) 0 ) 


The local truncation error in displaccinonls is obtained by inleg, rating 
the error in accelerations two times. A norm of enoi in di.splaienients 
([er? II can be expiessed as 


l|m-|| = (i + j^i) /.'IMI (2.''l) 

Since the above norm is an absolute quantity and problem itependent 
specifying an absolute value is not feavSible always. A relative I'rror 
norm ||enrei|| which is the ratio of error norm ||erj|| and noun of the 
displacements at the current step is used for step si/e eorieelion. /.c. 

f 





( 2 .‘) 2 ) 


Select the step size economically such that for early step, the loeal enor 
roughly equals the prescribed error tolerance. ‘iUH van he awhleved 
as follows: 


Define two parameters 7^ and ji such that 0 >, . j ami • I ami 


7lC < l|e}Tre/|| < 726 

Accept the step if (2.93) is satisfied. In practice the step swe adjustment 
is carried out if 
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1 Eiror is more than the upper limit of the toleiance ( 72 £) or 

2. Error falls below the lower limit of tolerance { 7 ic) ior fixed num- 
bci of times (say 5-10) consecutively This is to avoid costly 
factorization frequently needed whenever step change is made. 

New step size can be calculated as 

= (2.94) 

e-llrel 


2.6 Validation of the Model 

2.6,1 Introduction 

Before actually applying the model to three dimensional soil-structure 
iiiteraction problems which this investigation aimed to address, it is 
considered appropriate to validate the model with two-dimensional 
problems of similar kind. The set problems considered by Von Estorff, 
Pais and Kausel (1990) are chosen for validation purposes. The problem 
set consist of 1. A half space and 2. A half space with a vertical 
trench. To simulate transient loading a fast decaying Ricker's wavelet 
(Von Estorff, Pais and Kausel, 1990) has been chosen. The effect of 
model parameters such as number of Lanezos vectors (Section 2.4.5) 
considered, size of the mesh, element size (number of elements) have 
been studied and reported here. The implementation of the model for 
these problems has been discussed in the following sections. The effect 
of different (from that of the problem to be solved) spatial distributions 
of the load is also considered for the generation of Lanezos vectors, 
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It will be useful to undeistand the elfect o( spatial vaiiation ol tlu' 
load on the accuracy, and to a evolve a simpler eomputalion model 
for infinite domain independent of the loatlini^, 


2.6.2 FEM implementation 

A 2-dimensional finite element model with 4 noded iso|\uanu'lnc ele- 
ments is used in this implementation. A graded me.sh (['igiue 2.1)) with 
ratio of maximum element size minimum element .size (tenjdh/width} 
of 3,1 has been used. Since the problem is plane .sliain one, the 
elasticity matrix [E] becomes 


[E] = 


E 


(1 -I- i/)(i - 2u) 


1 - 1 / 1 / 

V 1--I/ 
0 0 


0 

0 

1 ■.>!' 
2 




and the strain displacement matrix (/,) is defmetl in (2.2.S). TIu’ .shape 
functions can be written as 


(1 - 0(1 - 0 



(1 + 0(1 ~v) 
(1 + 0(1 + 0 


i (1 - 0(1 + 0 J 


(S.ht)) 


Substituting these expression in (2.30) and integrating ov<«r the domain 
(area) of the element, stiffness matrix (A'"] can be written as 


[E\[D]dxdy 


(2.07) 
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wheie t IS the thickness of the element (taken to be unity) and A'' is 
the area of the element. In 2-dimensional domain each node will have 
2 degrees of freedom namely k, and v The matrix B in this case can 
be written as (it will be of oider of 3x8) 


[B] = [i] 


N 0 
0 N 


(2 98) 


In isoparametric formulation the physical domain is mapped onto para- 
metric domain using the same shape functions used for displacement 
approximation. The matrix operator [L] contains derivatives of depen- 
dent variables {u, v) with respect to physical coordinates (.r, y). Since 
the shape functions are expressed in natural coordinates C) 
derivatives with respect to physical coordinates have to be computed 
using Jacobian transformation from natural coordinate space to physical 
coordinate space. 

The derivatives with respect to natural coordinates can be expressed 
using chain rule as 


d j 6 dy 

d( ^dxdc dx ac 

d 4 - a dy 

dx di] ' dx dr] 


(2.99) 


The above equation can be expressed in matrix form as 


{ 


_d_ 


dx 

§}L 

OC 

a<: 

6x 



J 


dx 

dy 




_B_ 

dx 

dy 


( 2 . 100 ) 


where matrix |J] is called Jacobian of the transformation 
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Now the derivatives with lespecl to physical uKudin.Ue.s uin be ob 
tamed as 



The Jacobian in the above equation can be evaluated from (2.1(i) as 
follows 




(2.102) 


Where , indicates differentiation with icspect to variable follinvin}; it 
and Yj indicates summation over the number of uotle.s of the element, 
viz 






_ ^ 


.t'l + 




iUh.., 

d<, ‘b) 


L !i<Vi 
' il{ ‘'"I 


(2.103) 


Replacing the physical coordinates x and y by their natural coordinates, 
the stiffness matrix can be expressed as 


l^^'l = ' l'\of[E\[D]\\jy(;,li, (2, IIM) 

where 1|J1| denotes the determinant of Jacobian. 

Similarly the element mass matrix can be expre.ssed as 

M = /_|'/>|wriN|||y||<;c,;„ 


(j-iiiri) 
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The integration is carried out using Gaussian quadrature. The element 
load and damping matiices are not computed, rather global load vector 
in terms of global variables are computed since loading exists only at 
one node Since proportional (Rayleigh) damping is used the damping 
matrix is generated in the transformed (after applying the Lanezos 
transformation) domain. 

A value oi 9i =0 988 which produces small numerical damping is used 
in time integration The values of damping parameters a and 0 are 
taken to be 4% and 2% respectively. 


Ricker 's Wavelet 

{ 



Figure 2.3 Half Space Loaded by Vertical Rickers’s Wavelet 


2.6.3 Example 1. Elastic Half Space 

In this example an elastic half-space (see Von Estorff, Pais and Kausel, 
1990 ) subjected to a Ricker's wavelet (transient excitation) in vertical 
and horizontal direcfaons is investigated. The displacement histories 
at points 'A' and 'B' which are located at a distance of 20m, and 50m 
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Rickcj 's WjivclcJ 


> 



Figure 2.4: Half Space Loaded by llori/ontal Rickors's Wavcirt 
on the surface from the load aie monitored. 
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Figure 2,5, A Typical FE Mesh 


A schematic diagrams of the problem considered are shown in Fig- 
ures 2.3/ 2.4. Only one half of the domain is discretized for vertical 
loading and full domain is used for horizontal loading. A typical 
graded mesh with 20 number of elements in each direction shown in 
Figure 2.5, is utilized for vertical loading case. 
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Figure 2.6‘ Effect of Number of Vectors on Rel, Error Notm 


The parameters of the soil are as follows -- Younj^'s Motiulus ■ '.’.OO 
X IQ^kNlnx^, Poisson's ratio = 0.33, and density 2 The load 

function f{t) = (1 - where r - l/ir, 



Horzontal Displacement at A (m) Vertical Displacement at A (m) 
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Figure 2,1 • Effect of Number of Vectors on Vertical Displacement 



Figure 2.8, Effect of Number of Vectors on Horizontal Displacement 
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2.6.4 Number of Lanczos Vectors 

To study the effect of nimbcr of Lanczos vcciots considei lul foi com d i nalc 
transformation, the error norm in load represonlalion a.s dcliiu'd hy 
equation (2.78) is computed. Doth absolute and relative norms aie 
shown in Figure 2.6. The effect of number of Lane/o.s vecloj.s oji the 
accuracy of the results are shown in terms of vertical disiilaeeinent 
history at point 'A' in Figure 2.7. The effect can be .seen in the ease 
of horizontal displacement of half-space loaded by horizonlal Kiekt'r's 
wavelet in Figure 2.8. As expected error norm in lo.id lepussentalion 
reduces rapidly with increase in nuinbei of I.anc/.os vectms. .Snnilaily 
accuracy of the results improves rapidly as number of vectors Incroiisc’.s 
(reported in the following sections). It can be observed lhal 2(1 l.anc/os 
vectors could produce solutions of desired accuracy. 

Computational time 

To asses the computational resources required for Ihe lust two slap/'s 
the modeling scheme (Section 2.1), CPU limes are logged in and ana- 
lyzed in this section. The first two stages of the sdu’iiu’ (.Set (uui 2.1) 
namely the spatial discretization and I-anr/.os vector Iraiedonnaiion 
can be divided into the following steps - 1. IH* As.sembly, 2. l.ane/.os 
Vector generation, 3, Transformation matrix geiu'ration, d Writ' 
ing transformed system matrices and Lanczos vector.s (I/O) to disk. 
CPU times required for each of these steps obtained on l)i-t’ Al.1’1 1 A 
3000s600 are given in Ibblc 2.1. The s.inie liilarm.ilhm is slimvii in 
Figure 2.9. CPU lime required for slop I (PU assemlrly) rein.iins r-<m- 
stant with number of Lanczos vector bec.iusc lhi.s slop is In.lepniHlnnt 
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of Lanczos vector generation The CPU time lequired for each step 
per vector (siope of the lines in Figuie 2 9) remains same for steps 
2 and 4; whereas for step 3 which involves Ritz vector computation 
(Equation 2 81) increases. However, it can be observed that Total CPU 
time per vector (column 7) does not vary much 


No of 

FE 

Lane, vect 

Transform 

VO 

Total 

Total time 

Vectors 

Assembly 

generation 

matrix 



pei vector 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

5 

0.5232 

0 3867 

01376 

0.3670 

14145 

0.2829 

10 

0 5231 

0 6900 

0.3338 

0 7837 

2 3306 

0 2331 

20 

0.5241 

1 2171 

0 8042 

1.5460 

4.0914 

0 2046 

30 

05339 

1 8954 

1.8641 

23327 

6 6261 

0.2209 

50 

0 5348 

3.6246 

4 5193 

3,9138 

12 5925 

0.2519 

80 

0.5290 

75308 

10 2294 

6.3624 

24 6516 

0.3081 


Table 2.r CPU times for various stages of computation in sec 
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Figure 2.9; CPU time for various stages of coiniuaatum 


2.6.5 Size of the model 

To investigate the effect of the sixe of tho inoUol 
domain)/ four models of sizes 400x4()0m (small), 4l)(1x(4)()in (.'on.illl), 
600x60Qm (medium) and 800x80()iu (large) are couskleit'il iiv I'a'a’ of 
vertical loading. Number of elements in oai’h direction are taken to 
be 20. The time step is monitoicd carefully .so that Couutnt erilerlon 
(Higdon, 1992) is not violated. Lumped mas.'} mellunl is adojded 
for the analysis. Time integration is carried mil for 3see. Tlu' verti- 
cal displacement at point 'A' obtained with these models are 'ihown 
in Figure 2.11. The results obtained with tho models compare very 
well with those of Von Estoiff, Pais and Kausel (1990) obtained iisinp, 
boundary element methods. Only in small models, llie di-splaceinenls 
beyond 1.5 seconds do not match well, which coiiki be due lo the 
influence of reflected waves. 
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Similar study has been carried out with models 800x400m (small)^ 
1200x400m (smalll), 1200x600m (medium) and 1600x800m (large) for 
hoiizontal excitation. The number of elements used in hoiizontal and 
vertical directions are 40 and 20 respectively. Horizontal displacements 
at point 'K obtained with these three models are shown in Figure 2 13 
The results agree very well those reported by Von Estorff, Pais and 
Kausel (1990) (see Figure 2.12) except for those obtained with small 
meshes at times beyond 1 Ssecs 



Figure 2 10: Comparison of Vert. Displacements with Von Estorff's Values 




modeling dynamic SOIL-SIRUCIURI. INII.UAC MON 


^ 3e-06 

< 


C 

o 

e 

o 1e 06 


m 

(3 

e 

^ -1e-06 



I itloiH «<il 
Miiall 

jiiiuilM Mdiih 
Mod Mo'iti 
I ni (|0 Mo>ih 


0 0.5 1 1.5 !’ 

1 Imo ( 800 ) 

Figure 2 11: Effect of Mesh Size oti Vcrticn! DispLuoiiuMils 


E 

5e 06 

40-06 

< 

3e-06 

c 

0 ) 

P 

20-06 

Q> 

S 

1e-06 

s- 

Q 

0 

3 

C 

§ 

-1e-06 

o 

X 

-20-06 


-30-06 



Von I'olotd tiol. 

10 VOClOfil 
^^0 v(n.it)f;i 
00 VOl.llIf!) 


0.5 1 


Timo (aoc) 


Figure 2.12; Comparison of Flor. Displacements with Von Uunll's V.ilues 
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Figure 2.13 Effect of Mesh Size on Horizontal Displacements 


2.6.6 No of Elements 

In this section the effect of number of elements (element size) on the 
accuracy is reported. The computational domain of 800x800m size is 
discretized with varying no of elements from 20 to 40 in both the 
directions. The vertical displacement history at point 'A' subjected to 
a vertical load is obtained with these three mesh schemes and plotted 
in Figure 2.14. It can be observed that all the three models produced 
results of almost same accuracy and agree well with those reported by 
Von Estorff/ Pais and Kausel (1990). 
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Figure 2.14: Effect of Number of Elcmoiits 


dal Variation of the load 

n different loading schemes (shown in Figiu’e 2.15). an? 
rate the Lanczos vectors only namely 

at 10m right (node 22) of original po.silion (node 1) 
at 10m down (node 2) the original position (node 1) 
at 50m right (node 85) of original positiorr (node 1) 
at 50m down (node 4) the original po.sition (node 1) 
at 100m right (node 127) of original po.sition (node 1) 
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Figure 2 IS Spatial Variation of the load - Different Load Cases 


Large mesh with 20 elements in each direction is utilized for the 
study. Displacement histories at points 'A' and 'B' are obtained with 
20 and 40 Lanczos vectors. Vertical displacement history at point 'A' 
for vertical loading with 20 and 40 number of Lanczos vectors are 
shown in Figures 2.16, 2.17. It can be observed that results obtained 
with loading cases 2 and 4 match well; where as the ones obtained 
with other schemes do not match well with those obtained with actual 
loading scheme; It may be observed that if the loading scheme used for 
generation of Lanczos vectors is spatially farther (loading cases 3 and 
5) yields poor results and increase of vectors improves the accuracy. 
If the loading considered is close to the observation region, results 
obtained will be overestimated (loading case 1). 

It may be inferred that spatially different loading schemes can be used 
to generate Lanczos vectors and system matrices keeping in mind the 
following points. 


Vertical Displacement at A (m) 
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. The loading scheme should mil be voiy diK.'iviil (like UMdiiij., 
cases 3 and 5) from aclual loading, 

. The loading scheme should not close lo Ihe olwervaluin u'c.ion 
(like loading case 2) 
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Figure 2 16; Effect of Spatial Variation of the Loai) witli '.ill Voclors 
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Figure 2 17 Effect of Spatial Variation of the Load with 40 Vectors 


2.6.8 Example 2. Elastic Half Space with a Vertical 
Trench 

This example helps in understanding the behaviour and applicability 
of the scheme for non convex domains. A schematic diagram of half 
space with vertical trench adjacent to the load is shown in Figure 2 18. 
The geometry of the problem is same as in earlier example except that 
a vertical trench of 50m deep and 10m wide is located at 5m to right 
of the load. In this example also both vertical and horizontal loading 
schemes are considered. The displacement history at point 'K for both 
vertical and horizontal loading are shown along with those reported 
by Von Estorff, Pais and Kausel (1990) in Figures 2.19, 2.20. It can be 
observed that the results obtained by the present model agree well 
those obtained by Von Estorff, Pais and Kausel (1990), 
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RicKci ’s Wavt'lel 



Figure 2 . 18 : Half Space with trench Loaded hy Vortiuil RicIums's W.ivfh'l 



Figure 2,19. Effect of Trench on Vert. Displacemcnls due to Vert I o.ul 
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Figure 2 20‘ Effect of Trench on Hor Displacements due to Hor Load 

It would be important to mention that the solutions obey the causality 
conditions in all the situations investigated. This could be verified 
directly from the displacement histories reported. More over solutions 
obtained by Von Estorff, Pais and Kausel (1990) found to satisfy this 
requirement and the present solutions agree well with them 


2.7 Summary 

Based on the limited numerical investigation carried out the following 
conclusions can be drawn, 

1. A large mesh coupled with Lanczos vector transformation could 
represent semi infinite soil mass for computational purposes. 
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2. The present numerical scheme cai\ clfecUvely UamlU’ 
loading in time domain. 

3. The model can estimate errors due to aHn'din.Ue Iransloi inalion 
as well as temporal disci etizalion and control them. 

4. The error introduced due to coordinate transformation e.\[m- 
nentially decreases (see Figure 2.6) with increase in number ot 
Lanezos vectors. The computational cost per vectoi lemains al 
most constant and desired level ot accuracy can be obtained will\ 
sufficient number of Lanezos vectors. 

5. The system (transformed) matrices generated with a paitiudar 
spatial variation of the load can be used l<i analyze foi other 
loading conditions if the loading scheme.s are siinilai (see SeC' 
tion 2.6.7) 

6. A transformed model with few degree.^ of freedom (about 211 ■ '10) 
could simulate a physical model with large (aboiil H4tl) deg.rees 
of freedom. 

7. The present model can effocUvely analyze non convex iloin.iin.s 
also (vertical trench problems). 

8. The solutions obtained with present model obey the cau'nilily 
condition, 

9. With relatively small computational elfort sembinfinile half space 
problems with transient loading could be analysed in /im' t/omiiii. 

10. It would be important to emphasize (hat Ihe present ^.cheme is 
local both space and time. 



Chapter 3 


THREE DIMENSIONAL 
FOUNDATIONS 


3.1 Introduction 

Dynamic soil-structure interaction is concerned with study of sti'uc- 
tures founded on soils and subjected to dynamic loads directly applied 
to the structure or transmitted through the soil In its early stages, dy- 
namic soil-structure interaction dealt with external excitation problems 
like vibrations of machine foundations. More recently, important struc- 
tures, like nuclear reactor containment structures, tall buildings etc. 
excited through the motion of the soil in seismic areas, have received 
lot of attention of the researchers The flexibility of the foundation 
has a significant effect on the response of the structures and must be 
taken into account in the analysis. 

Soil structure interaction effects are studied usually in the frequency 
domain assuming hnear elastic behaviour of the soil. In this method 
the dynamic behaviour (responses) of the foundation (massless) along 
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with the soil is studied hrsl, Ihca the icspoasos .iloiip, Ihe louiul.uion 
structure interface are incorporated in the analysis ol the stun tine, 
Soil-structure interaction problems, wheie non lineai elhn h are iinpoi- 
tant, require a time domain analysis. Non-linear lonlait toiulilions, 
partial uplift, non-iineai material behavioui, and aibitrarily varying 
loads, are some of the usual silualions which iei|uire a diu'cl time 
domain analysis. 

Foundahons are usually massive, rectangular in shape and emlMulded 
in the soil. The usual concept of rigid fouiulation may not .ij'ply 
to the whole range of Ircqucncics of interest. C'on.seijoently the 
force displacement relationship of the foundtition, tiealinj; it to be 
rigid, through an impedance malrix may not de.scribe the behaviour 
accurately. 

The foundation soil system response depends on llie .shape and einbed- 
dment of the foundation, the rclalivo flexibility of soil iiivd foundation, 
spatial variation of the material propertie.s of the soil (mom hoinog.eMily 
and layering), and excitation. Because of the embeddmenl and timte 
rigidity of the foundation, in real life problems, the re.sptmses in the 
principal coordinated directions will be coupled. Any eompiehensive 
study of the foundation soil system should be able to ineorporale the 
above mentioned parameters in the model. A briel review of the litei- 
ature pertaining to the dynamic analysis of 3~dinu:nsioiuil foundations 
is presented in the next section. 
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3.2 Review of the Literature 

During the early part of this century machine foundations were de- 
signed by empirical models treating the system as single degiee of 
freedom systems (Balkan, 1969) Lamb (1904) initiated study of re- 
sponse of elastic half-space subjected to dynamic loads, Reissner (1936) 
obtained the response of harmonic disc load on surface He assumed 
a uniform distribution of stresses under the disc Quinlan (1953) re- 
ported solutions for varying contact pressure across the diameter He 
gave solutions for rigid base approximations. Sung (1953) developed 
solutions for circular contact area with three pressure distribuhons, 
namely, 1. uniform, 2. parabolic, and 3 pressure distribution cor- 
responding to rigid base. He gave solutions for different Poisson's 
ratios. These solutions can be used to compute the compliance func- 
tions of rectangular footing making use of the equivalent circular area. 
This equivalent area method provides solutions of reasonable accuracy 
for aspect ratio of two (E E. Richart, 1960). Sung (1953) developed 
mathematical expressions for rectangular footing resting on half-space. 
Kobori (1962), Thomson and Kobori (1962) obtained compliance func- 
tions for uniformly distributed load over rectangular surface area. 
Arnold, Bycroft and Warburton (1955) and Bycroft (1956) obtained the 
solution for rocking vibrations of a circular footing Actual stress dis- 
tribution beneath a footing does vary with operating frequency, and 
is not known apriori. A mixed boundary value problem, with pre- 
scribed displacements under a rigid footing and zero tractions over the 
remaining surface, was studied by Luco and Westmann (1971). Luco 
(1976) obtained impedance functions of rigid foundahon on layered 
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viscoelastic medium using integral equations approach. I’inile Ml'iiumU 
Method along with different absorbing boundaiics wore employed by 
several investigators to obtain the impedance functions o( suilace and 
embedded foundations. Lysmer and Kiihlemcycr (1969) developetl 
an absorbing boundary consisting of discrete viscous tlampeis and 
Rayleigh wave dampers. They obtained compliance fiincUons for 
strip and axisymmetne foundations. Dasgupla and Kainesvvaia Ran 
(1978) extended the methodology to study plane, axisymmelric and 
three dimensional problems. They considered non-linoar maleiial be- 
haviour and non-homogeneity. Lysmer and Waas (1972) (.leveloin'd 
an absorbing boundary which transmits Love and Rayleigh waves. 
This boundary was applicable for layered media with hori/ontal layers 
extending to infinity in lateral directions resting on hard strain in In 
frequency domain. Katisel (1974) extended the same for axisymmelric 
problems and obtained impedance functions for different modes. 

The Boundary Element Method (BEM) which automatically includes 
infinite region m the modeling have been extensively used to study 
the soil-structure interaction problems. The frequency tiomain for- 
mulation was first used by Dominguez (1978) to obtain impedance 
functions of rectangular foundations testing on or embedded in a 
viscoelastic half-space. Impedance functions of foundations resting 
on non-homogeneous and layered soils were obtained by Abascal and 
Dominguez (1986). Ahmad and Banerjee (1988) employed quadratic el- 
ements in frequency domain to obtain compliance functions. Karabalis 
and Beskos (1984) and Karabalis and Beskos (1986) used a time -domain 
BEM and obtained impedance funclions for rigid surface and embed- 
ded foundations. The response of flexible foundations was invesligaled 
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making use of domain methods like Finite Element (FE) and Finite Dif- 
ference (FD) foi foundation and Boundary Element (BE) for soil Iguchi 
and Luco (1981) obtained the response of lectangular foundations on 
elastic half-space. Karabalis and Beskos (1985) and Gaitanaros and 
Karabalis (1988) obtained response of surface and embedded flexible 
foundations using time domain boundary element method and FEM 
for the foundation, Gucunski and Peek (1993) piesented response of 
circular foundations using soil stiffness matrices for layered soils eval- 
uated through Green's functions and stiffness matrices of plate using 
finite difference energy method. 

Gazetas and Tassoulas (1987b), Gazetas and Tassoulas (1987a) devel- 
oped approximate models for stiffness and damping parameters for 
horizontal and rocking vibrations based on results obtained through 
numerical computation. Wolf (1988) developed approximate models by 
adding an additional degree of freedom. He obtained stiffness, damp- 
ing and mass parameters by curve fitting techniques to the results 
obtained through exact numerical computation. 

The existing methodologies to solve the 3dimensional foundation soil 
systems fall in one of the categories. 

1, Finite element discretization with energy absorbing boundaries 
to model the far field. 

2, Boundary element method in frequency domain ignoring the 
foundation flexibility. 

3. Boundary element method in time domain- 

4. Hybrid methods - FE discretization for foundation and BE for 
the rest of the domain. 
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Category 1 methods account the energy radiation approxiiiiiitcly. Most 
of these models, like those of Kausel (1974) and Lysincr and Kulilemcyer 
(1969); can be used in frequency domain only. Category 2, .1 airtt 4 
methods can model the infinite extent exactly. The lime domain 
variants of these BE methods, being global in space and time, need 
enormous computational resources. They can not handle flexible 
foundations directly. Category 4 methods also suffer the tlisadvanlagi'.s 
of the BEMs in time domain. 

The present investigation uses a finite clement model for spatial dis- 
cretization (both near and far field) circumventing the need for an 
energy absorbing boundary. The response of the sy.stein is obtained in 
time domain by direct time integration after transforming the equations 
from onginal space to Lanezos space. The investigation thus uses in 
integrated approach to solve coupled motion wherein all niode.s could 
be present. 


33 Application of the Finite Element Method 

Analytical solutions to dynamic soil-structure inlcraelion problems can 
be obtained for situations with idealized geometry of the foimdalion, 
linear material behaviour and periodic loading. Numerical iec'hnii(Ui'S, 
with advent of high speed computers, provide methodology to solve 
complex field problems. The Finite Element Method (FEM) and Bound- 
ary Element Method (BEM) are some of the powerful methods which 
employ discretization of the domain. The FEM, as pointed out in 
Chapter 1, is ideally suited for problems with flexible foundation sys 
terns, non-homogeneous and nonlinear material behaviour and for 
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time domain analysis, if the effect of far held is incorporated properly. 
In this investigation the spatial disci ctization is cairied out using FEM 
considciing a Imgc domain of the half-space, Lanezos transformation 
techniques (described in Section 2 4) are applied to system of equations 
of motion resulting from the spatial discretization The transformed 
equations are solved in time domain with an adaptive diiect integiation 
method described in Section 2.5 3. 

3.3.1 Finite Element Model 

An eight noded hexahedron isoparametric element is consideied as 
basic element. Considering the range of foundation thickness and dis- 
placement continuity with underlying soil, it has also been discretized 
with same eight noded elastic elements. The governing equation of 
motion ((2.49)) in which the system matrices (A/j, [CJ and [A] are as- 
’ sembled from the individual element matrices. The elements mass, 
damping and stiffness matrices are obtained as described in the fol- 
lowing section. 

Element stififness matrix 

An eight noded isoparametric element is shown in global and local 
coordinate systems in Figure 3.1. Global node numbers of a typical 
element are shown in italics and local node numbers for isoparametric 
element are shown in roman numerals. Since the local coordinates 
range from -1 to +1, the element is a cube with a side of 2 units long 
and origin of the system at its centroid. The actual element could be 
of any shape with six-faces. The coordinate mapping which relates 
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the two coordinate system can be expressed as (Equalion 3.2) 



T 


(a) Global coordinate sysicm (b) I^Kal cmndinule system 


Figure 3.1 8 Noded Hexahedron Element in Global and Local Coordinate System 


X=:[N.1{X.} 


(.3 1) 


where x and x, are the vectors of independent vaiiahle.s, and cotudh 
nates for the node i respectively; (N.j is a diagonal .sub inaliix of oidor 
equal to the number of space dimensions (3). Each of Its entry otpialH 
to mapping function of the node. The subscript i vark'.s from 1 to H. 
In (3.1), 


X = [a; 1/ 2 l'^ 


[N.l = 


iV, 0 0 

0 /V, 0 

0 0 yv, 


X| — X‘( l/i Zf 


( 3 . 2 ) 
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The mapping fiinchons can be written as 

A'. = (l + CiC)0 + J/,7,)(l+eO (3 3) 

where C, V>^ C< ate coordinates of the node ? and take eillier +1 
or -1 

The displacement vector u can be expressed in terms of its nodal values 
using the shape function expressed in the local coordinate system of 
the element Since isoparametric formulation is being adapted the 
shape functions are same as the coordinate mapping functions of the 
element. Hence, 


u = [Ni]{u,} 


(3.4) 


The differential operator defining the stress strain relations can be 
expressed as 


\L] = 


dx 

0 
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dx 
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dx 
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dz 
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dx: ^ 


(3.5) 


The strain nodal displacement matrix B is given by 

\B] = [im 


(3.6) 


The differential operator [L\ contains derivatives with respect to phys- 
ical coordinates. The derivatives have to be computed using Jacobian 
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transformation from natural coordinate space to physical coordinate 
space. The Jacobian can be expressed as (detailed description is jpven 
in Section 2.6.3), 



I]!A^i,C‘^'t 


(3.7) 


The ' indicates differentiation with respect to the variable following; 
the symbol and the symbol ^ indicates summation ovei the mimber 
of nodes of the element (8 here), viz. 


Yj + ^a-'a + -I- 



(3.8) 


The elasticity matrix [/5J for three dimensionnl solid can be e.xpie.s.sed 
as, 


B] 
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The stiffness matrix of the element now can be expre.s.sed 

= £ j'' jjmmymi'ii 


(••i.'i) 


(3. Ill) 


where IIJ^ denotes the determinant of Jacobian. The slifiness matrix 
is evaluated using a (2x2x2) Gaussian quadrature. 
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Element mass matrix 

Element mass maliix can be expiessed as 

yy r PIN]’ |w||| (3.11) 

As already mentioned, systematic lumping method is adopted for mass 
matrix Only the diagonal entries of the element mass matrix are 
evaluated using (2x2x2) quadrature and they are scaled so that the 
element mass is conserved. 

Damping matrix 

Element damping matrices are not evaluated since Rayleigh damping 
(proportional) is used, the damping matrices (global) in the transformed 
domain are computed directly. 

Force vector 

Since very few elements are loaded, element load vectors are not 
computed, rather the global load vector is computed. Most of the 
computations are done using point loads. 

Linearly varying Young's modulus 

Soil medium with hnearly varying Young's modulus has also been 
considered for the analysis. The stiffness matrices of the elements are 
computed using an effective Young's modulus. Young's modulus in 
each element is assumed constant and its magnitude is taken to be 
its value at the centroid of the element. The variation of Young's 
modulus with depth (y coordinate) may be expressed as 
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where Et, Ey are the modulii at the surface ( /.f. y - 0) and at depth i/; 
a is a constarrt which equals to the slope of the vaiiation. Respnuses 
are obtained for different values of a. 

Layered soil 

Horizontally layered soils have been considered for the analysis. 
Young's modulus of each layer can have different value.s. 'I'lie dis- 
cretization is done such that the element boundaiies lie alonj? the 
interface between the layers. The stiffness matrices of the elemenl.s 
are computed using the Young's modulus of the layer in which they 
he. A two layer system is investigated with different ratios of Young's 
modulii and different top layer thickness. The ratio of the thickness 
of the top layer to halfwidth of footing and ratio of Young's inothilus 
of bottom soil to that of top soil are used to present the res[)onhes. 
Poisson's ratio and density of the layers arc kept same for all layers as 
they do not significantly effect the responscDasgupta and K.uneswara 
Rao (1978). 

Rigid Bedrock at Finite Depth 

Soil of finite thickness lying on a rigid bedrock is also iuve.s(i}'ated. 
The displacements at bottom of the soil stratum are restrained in this 
analysis. A non-dimensional thickness parameter, which is defined 
as the raho of the thickness of the soil layer lo the halfwidtli of the 
footing is used to present the responses. 
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Footing cmbcddmcnt 

In gcncial footings arc embedded in the soil To study the effect 
footings at diffeicnl embeddnients are consideied, Loading is applied 
on the top surface of the foundation. A non-dimeiisional ratio of 
depth of embeddment to halfwidth of footing is used to piesent the 
responses 

Foundation flexibility 

Foundation flexibility has significant effect on responses. A non- 
dimensional flexibility ratio (reciprocal of !]gidity)is used to chaiacterize 
the relative flexibility of the foundation. Vaiying the Young's modulus 
of the foundation material, different flexibility ratios are obtained. 
The effect of foundation flexibility in terms of the variation of the 
displacements along width and/or length of the foundation is used to 
present the results. The effect of embeddment of the foundation on 
flexibility is also investigated. 


3.4 Presentation of Results 

Amplitudes for vertical, horizontal and rocking (rotation about z) modes 
of vibration are obtained by applying harmonic loads/moments in 
each of the direction independently. Integration is carried out till 
the transients die and steady state is reached. It is found that if the 
integration is carried out for about 8 to 9 times the periods of the load, 
steady state could be achieved. The amplitudes are obtained from 
the displacement histories. Compliance functions for rigid foundation 
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are computed for vertical, coupled hodzonUil anti lotalinnal iiukIl's 
of vibration for frequency latios ranging from 0 D It' '-I 0. A non 
dimensional compliance function is defined for each nuulo ol the 
vibrahon for any given frequency ratio flo as, 


^ _G,L^ 

Cii p ^ 




^M>yy- u 


i>y 


and 


uB 



(:ii:i) 


{^.U) 


where C^y, and are the compliance functions for ver- 

tical, horizontal along x direction, rotation about j axis anti loisional 
directions respectively; v is the vertical amplitude due to vertical load 
Py, u is the horizontal displacement along x axis due the loatl acting 
along a; direction, is the rotation about z Axis due to tlu' momenl 
about 2 axis, ipy is the torsional displacomenl (about // axis) due to 
the torsional moment My about y axis, tu is the frequency of excitation. 
Torsional compliance functions arc obtained only for the eonqiarison 
purposes 



where p, is the mass density of the soil. 
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Complinnce functions for coupled honzontal and rocking modes are 
also ohlajned for rigid foundation 




C'srt/r 




GJ.D , 

y>. 

G,LB 


G^iZT Gx^z 


(3.16) 


Here the compliance functions C^,zr is the compliance relating the 
lolalion iJk about ^ axis due to Horizontal load Px along x axiS/ 
similarly C^^z relates the the horizontal displacement v along i axis 
due to moment Mz about x: axis. 


3.4.1 Validation of the model 

A typical 3-dimensional foundation excited by harmonic loads and 
moments in vertical, and coupled horizontal and rotational modes 
independently, is analyzed with the present model. Problem definition 
and a part of the discretized domain used for vertical mode are shown 
in Figures 3.2, 3,3 respectively. Node numbers are shown in roman 
numerals and element numbers are in italics. L, B and d are length, 
width and thickness of the foundation respectively; and h is the depth 
of embeddment. Due to symmetry one quarter of the model for 
vertical mode and one half of the model for horizontal mode is chosen 
for discretization. The geometry and the material properties of the 
problems analyzed are as follows — 


Vertical Mode: 
Size of the model 


30x30x30m 
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Number of Elements ^ 

Young's Modulus of Soil E, = 
Poisson's Ratio Soil Ua ~ 

Unit weight of Soil pa = 

Half Length of Footing L = 
Half Width of Footing B 
Thickness of Footing d = 

Depth of Embeddmenl h = 


12x10x12 

47.09MPa 

0.33 

23.52 kN /m'’ 
2.0m 
2 Om 
1.0m 
1.0m 


Horizontal Mode: 

Size of the model = 60x30x30m 

Number of Elements = 24x10x12 

Young' s Modulus of Soil E^ = 47.09MPa 

Poisson's Ratio Soil t/, = 0.33 

Unit weight of Soil Ps = 23.52 kN /m ’ 

Half Length of Footing L = 2.0m 

Half Width of Footing B = 2.0m 

Thickness of Footing d = 1.0m 

Depth of Embeddment h = i.Om 


Compliance functions obtained with present model are compared vvilli 
those reported by Mita and Luco (1989) and arc shown in Ligiires .14, 
3.5, 3.6 and 3.7. The results obtained compare well with those reported 
by Mita and Luco (1989) 



3 4 Presentation of Results 


95 



Figure 3 2 Problem Definition for 3dimensional Embedded foundation 



Figure 3.3: Finite Element Discretization forSdtmensional Embedded foundation 
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Figure 3.4‘ Comparison of Vertical Compliance Function Values with Mila's Val- 
ues 



Figure 3.5: Comparison of Horizontal Compliance Function with Mita's Values. 



Vyy 
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Figure 3.6, Comparison of Rotational Compliance Function with Mita's Values 



Figure 3.7: Comparison of Torsional Compliance Function with Mita’s Values, 
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3.4.2 EJtfect of Embeddinent 

Compliance for vertical, horizontal and rockinj; inoilos for vaiious 
embeddment ratios ranging from 0.5 to 1.5 (ratio ol doplh ol foiiiulalion 
to halfwidth of footing lifB) are computed and pieseiitod in t-ip,nu‘-, 
3.8 to 3.16. It can be observed that with inciease in embeddinent ratio, 
the compliance functions decrease throughout tlie frequency lange. 
However the reduction in the compliance functions vanes with mode 
and aspect ratio of the foundation. The reduction in Ihv' valnc.s nl 
compliance functions for vertical modes with embedilment is small (nl 
the order of 10%) except for static case (TIguies 0.8 to O.ID). In I'lisc of 
horizontal mode the reduction is of the order of 50% (I'ip, tin's O.l 1 to 
3.13). Similar order of reduction can be observed for rotational morlc 
of vibration (Figures 3.14 to 3.16). 

For vertical mode of vibration the reduclion in compliance flint lion 
with increase in embeddment ratio is more (or square fount lation 
than for rectangular one (Figuies 3.8 to 3.10). Similar behavitme Is 
observed for horizontal mode of vibration (Figures 3.11 lo 3.13), Ftiv 
rotational mode of vibration reduction in compliance funelion i.s intire 
for rectangular foundation than square one (Figure.s 3,M It) .l.lh). 
Embeddment of the foundation does not significantly alter Ihv' resonant 
frequencies for any of the modes of vibration. 
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Figure 3 10. Effect of Embeddment on Vertical Compliance TmiCtion for Rect* 
angular Foundation (L/B = 3) 
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Figure 3.11 Effect of Embeddment on Horizontal Compliaiuc fTinclion for 
Square Foundation 
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Figure 3 12 Effect of Embeddment on Horizontal Compliance Function for Rect- 
angular Foundation (L/B = 2) 



Figure 3,13’ Effect of Embeddment on Horizontal Compliance Function for Rect- 
angular Foundation (L/B = 3) 
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Figure 3.14 Effect of Embeddment on Rotational Cornpll.'iiKO (or 

Square Foundation 



t 


Ut\ 

Km 

x 




i I 

1.G 1.0 


Figure 3.15 Effect of Embeddment on Rotational Compliance l“i,nction for Rod. 
angular Foundation (L/B -- 2) 
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Figure 3.16 Effect of Embeddment on Rotational Compliance Function for Rect- 
angular Foundation (L/B = 3) 

3.4.3 Effect of Aspect Ratio (L/B) 

The Effect of Aspect ratio of the foundation on the compliance function 
at different embeddment rahos are investigated. Compliance functions 
for vertical, horizontal and rocking modes are shown in Figures 3.17 
to 3.24. 

It can be noted from Figures 3.17 to 3.24 that aspect ratio does not alter 
the resonant frequencies in all three modes of vibration. For vertical 
mode of vibration, compliance function (Figures 3.17 to 3.19) increases 
with aspect ratio by more than 100% from square foundation to rect- 
angular foundation with LjB = 3 at resonance. At high frequencies 
the effect of aspect ratio is small. For horizontal mode of vibration 
compliance function (Figures 3.20 to 3.22) increase with aspect ratio 
by about 90% from square foundation to rectangular foundation with 




104 


THR EE DIMLNSIONAI f OUNDAI IONS 

LjB = 2 at resonance. In case of rolalional vibration llu’ a.spiHl lalio 
does not have significant effect (Figures 3.23 to 3.24 a.s in lase llu* ol 
other modes. Maximum increase in the compliance (unction is ol the 
order of 20%. 



Figure 3 17: Effect of Aspect Ratio on Vertical Compliance I unction for 1 mbecl- 
ded Foundation (h/B = 0 5) 
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Figure 3 18 Effect of Aspect Ratio on Vertical Compliance Function for Embed 
ded Foundation (h/B = 10) 
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Figure 3. 19’ Effect of Aspect Ratio on Vertical Compliance Function for Embed 
ded Foundation (h/B = 1.5) 
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Figure 3 20 Effect of Aspect Ratio on Horizontal Compliance l uiu lt<in lor l-m 
bedded Foundation (h/b = 0.5) 
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Figure 3 21 Effect of Aspect Ratio on Horizontal Compliance (•unction for Ctn 
bedded Foundation (h/B = 1.0) 
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Figure 3.24: Effect of Aspect Ratio on Rotational Compliance f-unction lor I nv 
bedded Foundation (h/B ~ 1.5) 


3.4.4 Effect of Foundation Flexibility 

A non dimensional relative flexibility is delincd to cliaracli'ii/o tlu' 
foundation flexibility as (Gucunski and Peek (19^X1)), 


where Ep, are Young's modulus of the foundation maloiial and soil 
respectively. 

Displacement amplitudes are computed for .square I'lnlH'ddod foimda’ 
tion with flexibility ratios Fv of 0 125, 0.0125, and 0 0(1125. Normali/t'd 
amplitude is defined as the ratio of the ampUuidvs at any point on 
the surface of the foundation to that at the center. A nornitili/eil 
amplitude of 1.0 all along the surface indicate.s a rip, id I'ouiulalion. 
The absolute and normalized variation of the ampliliules along the 
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width of the foundnlion at diffcicnt frequency ratios for vcihcal and 
lioii/onlal modes of vibiation are obtained 

The vaiiation of absolute amplitudes for vertical mode along the nor- 
malized width (latio of distance t of the location of the point from 
the center to halfwidth of the foundation) is shown in Figure 3.25 
Absolute amplitudes rapidly deciease along the width for high flexi- 
bility ratio, wheieas they do not change much for low flexibility ratios 
Howevei beyond the edge of the foundation flexibility does not have 
significant effect. Normalized amplitudes at diffeient frequency ratios 
and flexibility latios are shown in Figure 3.26. It can be observed 
that increase m frequency ratio has little effect on the normalized 
displacements Normalized displacement variation with different em- 
beddment ratios are shown in Figme 3.27 Increase m embeddment 
ratio increases the effective rigidity of the foundation, t.e. for the same 
flexibility ratio increase in the embeddment reduces the variation of 
the normalized displacement. 

Variation of absolute and normalized horizontal amplitudes for differ- 
ent flexibility ratios are shown in Figures 3.28 and 3.29. It can be seen 
that flexibility ratio does not have significant effect in horizontal mode 
of vibration. 




x/B 


Figure 3 25- Effect of Foundation Flexibility on Absolute Vcrtk’.il Anipliliidc 
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Figure 3 26 Foundation Flexibility vs Normalized Vertical Ainplilii<lo 
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Figure 3.27 Embeddment vs Normalized Vertical Amplitude 



Figure 3.28. Effect of Foundation Flexibility on Absolute Horizontal Amplitude 




112 


THREE DIMENSIONAl [OUNOAIIONS 



Figure 3,29 Effect of Foundation Flexibility on Normaliecd I loii/ontjl Anipliluiln 


3.4.5 Effect ot layered Soil 

Compliance functions of square embetWed fouml.dion., ivslin,; 

two-layered soil system are computed, liffeci of lop |.,yo,. tluckncss 

and ratio of Young's moduli, of the two layers are invesligated. No,,- 

dimensional thickness parameler '/■, and modular lalio A/„ are ilelined 
as 


rp = 


M, 


B 

St 


p.lH) 




where i is the thickness of top Idyer- F mrl /<' i v 
of ton larro c4 1. y • ‘t- /rfc rU’c Ilu' Yoimp's niodulii 

top layer and bottom layer respectively. 

particular modular ratio are shown in M«nres 3.30 
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52. Il can be seen that the effect of layering is pronounced ( 
?diiclion in compliance function foi M, >10 and inciease for 
1.0) foi small values of thickness latio It can be said that small 
ness of the top layei effects the responses consideiably In the 
2 mentioned hgures the 'half-space' solution lepiesents response 
homogeneous half-space with modulus equal to It can be 
ved that beyond a thickness ratio of 10 layering does not affect 
2 sponse i e. top soil whose thickness equal to 10 times halfwidth 
)ting will only influence the lesponse. The presence of a soft soil 
(M, < 1) under-lain reduces the resonant frequencies, whereas a 
oil > 1) inciease the resonant fiequencies, A layered soil sys- 
vith higher modular ratio with small thickness of the top soil do 
xhibit resonant frequencies. Figures 3.34 and 3.33 show the effect 
idular ratio at different thickness ratios. It can be observed that 
of modular ratio is strong at low frequencies and low thickness 




Froquor^cy DqIIo 


Figure 3.30; Thickness Ratio vs Vertical Compliance ( unction (M, 0 !i) 



Figure 3 31: Thickness Ratio ys Vertical Compliance Function (M, 2 0) 
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Figure 3 34 Modular Ratio vs Vertical Compliance I uiu lion ( I’, H 0) 

Horizontal compliance functions for different thickness ratios al mod 
ular ratios of 0,5, 2.0 and 4.0 are presented in Inj’ures to It IV/. in 
horizontal mode also a small thickness of the lop layer ellects llie ie» 
sponses strongly. For low T,., resonant frequency of Ihe .system u'lUices 
for Mr < 1 and increases for M, > i. Beyond a thickness ol ratio ol H, 
layering does not have significant effect on the coiniiliance functions. 
Compliance functions for different modular ratios are .slunvn in I’iji;’ 
ures 3.38 and 3.39 for thickness ratios of '2 and H lespeclively, It lan 
be seen that unlike vertical mode the modular ralio does not havi- any 
significant effect beyond frequency ratios of 0.0 even for tow thickness 
ratios. 

Effect of layering on rotational mode of vibration is presented in 
Figures 3,40 to 3,41, As in the case of horizontal modi’, rolalional 
compliance functions are significantly influenced by small thick lu-ss of 
top layer. However the reduction for systems wllli Mr > 1 or increa.se 
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for systems with 71/, < l, in the compliance functions is small (of the 
order of 30%) 



Frequency Ratio 


Figure 3 35 Thickness Ratio vs Horizontal Compliance Function (Mr = 0 5) 



Figure 3.36; Thickness Ratio vs Horizontal Compliance Function (Mr = 2,0) 
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Figure 3 37 Thickness Ratio vs Horizontal Compliance Fiim tion (M, A 0) 
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Figure 3.38. Modular Ratio vs Horizontal CompHance l untliori ( )', 2 0) 
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Figure 3.41 Thickness Ratio e/s Rotational Compli<l(tce Fiioclion (M, 2.0) 


3.4.6 Effect of Linearly Varying Young's Modulus 

Effect of linearly varying Young's modulus on compliance functions of 
vertical and horizontal modes of vibration is reported in Ihis section. 
Vertical and Horizontal compliance functions aie shown lor diffi’reiU 
values of the parameter a in Figures 3.42 and 3.43. U can observed 
that linearly increasing modulus reducc.s the compliant' /iiiu liou.s. 
The reduction is more pronounced in horizontal mode than In veilical 
mode of vibration. In horizontal mode of vibration an iuciease in the 
value a increases the resonant frequency. 
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Figure 3 42' Effect of Linearly Varying Modulus on Vertical Compliance Function 



Frequency Ratio 

Figure 3.43 Effect of Linearly Varying Modulus on Florizontal Compliance Func- 
tion 
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3.4.7 Effect of Rigid Bedrock at Finite Depth 

Vertical and Horizontal compliance funclions for finite thick (I)) im- 
pressible layer resting on rigid bedrock are prcsenled in this section 
For different values of the ratio DJ13, where 1) is the thickness of the 
compressible layer, compliance functions for vertical and hori/.ontal 
modes of vibration are shown in Figures 3.44 and 3.45. It can be 
observed that vertical responses are strongly inlluenced, wherea.s the 
horizontal responses are less influenced by the presence ol liarcl roik 
at finite depth. 



Figure 3 44: Effect of Rigid Bedrock at Finite Depth on Vortical Conipli.inr o 
Function 
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Figure 3 45 Effect of Rigid Bedrock at Finite Depth on Horizontal Compliance 
Function 

3.4.8 Summary 

A summary of the conclusions based on the results presented in the 
earlier sections are presented in this section. 

1. The solutions obtained by the present model match well with 
those obtained by the time domain BEM method. 

2. Reduction in compliance functions due embeddment is more in 
horizontal and rotational modes than in vertical mode. Embed- 
dment does not significantly alter the resonant frequencies of the 
system. 

3. Increase in aspect ratio (IVB) causes an increase in the compliance 
functions in all modes. The effect of aspect ratio is small in 
rotational mode. However at higher frequencies the effect of 
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aspect ratio is small for all modes of vibralion. 

4. Foundations behave less rigid (can be quanlifiod in ienns of iho 
variation of displacement across) in verlicnl niocio than in hoi izon- 
tal mode for given geometry and mateiial properties (floxibllily 
ratio). Increase in the depth of embeddment increases Iho rigullly 
of the foundation (in vertical mode) for a given fioxibilUy ratio. 

5. Compliance functions are strongly influenced by small thickness 
of the top layer. Even at low frequencies/ soil layer below 10 litiu’s 
the half width of the footing docs not influence Ihe behaviour. 
Due to layering resonant frequencies changes. Soft boltom layiiis 
reduce the resonant frequencies whereas the hard bottom layer 
increase the resonant frequencies of the syaletn. 

6. Linearly increasing modulus reduces the compliance functions. 
The effect is more pronounced for horizontal mode than vortical 
mode. Linear variation of the moduhi.s incieases the resonant 
frequencies for horizontal mode, whereas vertical mode doe.s not 
exhibit any resonance. 

7. presence of hard rock is strongly fell in vertical mode than 
horizontal mode of vibration. 



Chapter 4 

FILE FOUNDATION SOIL SYSTEMS 


Pile foundation aie usually adopted in situations where — 1. The 
lop soil is highly compressible, 2. It is not strong enough to take 
the loads through block foundation, 3. It is prone to liquefaction 
4. To transfer large magnitudes of lateial loads to the soil. The 
response of single pile and pile groups subjected to dynamic loads 
have received a lot of attention over last couple of decades. The 
development of analytical solutions is accelerated due to growing use 
of pile foundations in traditional areas such as buildings, machine 
foundations, nuclear power plants and offshore structures. 


4.1 Introduction 

The response of piles and their load carrying capacity depends on the 
pile soil interaction. This interaction modifies the stiffness of the pile 
and generates damping through energy radiation (geometric damping) 
due to semi-infinite geometry of the soil and energy dissipation due 
to material damping. Often the piles are used in a group, which 
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introduces pile-soil-pile internclion. The coinhiiu'd cUcct ol pile <'ap 
and pile group determines the response of the pile loiindalion syMein 
Due to the complex nature of soil-structure (pile and/or cap) iiUciaciion 
(see Chapter 3) analytical solutions are difficult to obtain 


4.2 Brief Literature Review 

Novak (1974) was the first to attempt to analy/.e tlu- dynamic rc 
sponse of single pile using continuum theory, Ilo moilclcd Ihc pilt* 
as an assembly of beam elements and soil roacliouH on i( vvert' oli- 
tamed assuming a plane strain condition (infinitely long pile ). T'lie 
same approach was extended by Novak (1977) for llu* analysis of pile 
embedded in layered media. Nogaml and Novak (1976) and Novak 
and Nogami (1977) developed closed form solulions/ for verlical aiul 
horizontal modes of vibration of an end bearing pile, making use of 
rigorous three dimensional theory. Novak and Ifowell (197H) obliiineil 
the responses of the pile m torsional mode of vibration. Novak and 
Sheta (1980) have investigated the response of piles including, the con- 
tact effects using a mass less layer with reduced modulus .noiuul Ihe 
piles. Novak and Sheta (1982) obtained ihc pile group responses using 
interaction factors. 

Blaney, Kausel and Roesett (1976) studied pile vibration using, (mile 
element method with consistent energy absoibing boundary. Kuhle- 
meyer (1979b), and Kuhiemeyer (1979a) studied vertical aiul hon/on- 
tally loaded piles by finite element method, ffo developed a method 
for axisymmetrical solids loaded non-axlsymmclrically, making use of 
finite elements and energy absorbing boundary. Dobry, Vincenli*, e/ u/. 
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(19H2) studied the horizontally loaded piles making use of finite ele- 
ment codes available and piesented some approximations for stiffness 
and damping ratios. Wolf and von Arx (1978) extended the finite el- 
ement analysis to pile groups. Waas and Haitmann (1981) developed 
a methodology to analyze pile groups m which pile stiffnesses were 
obtained using a Fourier expansion for axisymmetrical pile groups and 
finite elements along with consistent energy absorbing boundary for 
soil This was extended by Tyson and Kausel (1983) for pile groups 
with embedded pile laft. Kaynia (1982) developed a three dimensional 
continuum model for pile groups and obtained dynamic interaction 
factors and group efficiency latios 

Sen, Kausel and Baneijee (1985) developed boundary element formu- 
lation to obtain the steady slate response of piles and pile groups 
embedded in homogeneous and non-homogeneous soils. Nogami 
and Kongai (1986) and Nogami and Kongai (1988) developed discrete 
models for pile soil interaction for vertical and horizontal modes and 
obtained responses in time domain. Mamoon and Banerjee (1992) 
developed a direct time-domain and transformed time domain models 
for dynamic pile soil interaction. They used step-by-step integration 
scheme using approximate half-space half-space integral formulation 
to obtain the soil reactions. Cheung, Tham and Lie (1995) used direct 
boundary element formulation to obtain the response of horizontally 
loaded pile in time domain. 
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4.3 Present Model 

The procedures for frequency domain analysis of pile foundations 
are well investigated Extensive literature is available on dynamic 
response of single piles. However, usually, piles are used in gioups 
with a cap. As mentioned in the review, the use of dynamic interaction 
factors to obtain the group response is atmost approximate in nature. 
It IS rather prudent to visualize the pile group with cap as a raft 
foundation supported by piles. Here in this investigation the pile 
group has been analyzed along with pile cap (or as block foundation 
supported by group of piles) in an integrated manner, rather than 
the substructure approach (treating soil and pile separate) adopted by 
earlier investigators. To appreciate the effect of different parameters 
on the response of piles a bnef study has been presented for single 
pile response. The computational model used for this investigation is 
essentially the same utilized for block foundation analysis described 
in detail in Chapter 3 (Section 3). The piles are modeled using three 
dimensional 8 noded continuum elements. The time domain solutions 
obtained by the present model for single pile have been validated with 
existing solutions. 

4.4 Validation of the Model 

A typical pile vibrating in vertical and lateral modes due a transient 
excitation has been chosen for validation. The problem definition is 
given in Figure 4.1. The excitation is a triangular pulse and is shown 
in Figure 4,2, where P^ax is maximum force and duration of the 
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oxciUilion is A/. /, d are the length and side of a square pile which is 
a good nppioximnlion for a circular pile The geometry and material 
piopeilies of the problems considered in the analysis are as follows 


Veitical Mode: 

Size of the model 
Number of Elements 
Young's Modulus of Pile Material Ep 
Poisson's Ratio Pile Material 
Density of Pile Material pp 
Young's Modulus of Soil E, 

Poisson's Ratio Soil 
Density of Soil 
Length of Pile I 
Diameter of Pile d 

p 

J max 

Horizontal Mode: 

Size of the model 
Number of Elements 
Young's Modulus of Pile Material Ep 
Poisson's Ratio Pile Material Up 
Unit weight of Pile Material Pp 
Young's Modulus of Soil Es 


= lOOxlOOxlOOm 

= 15x15x15 (in X, y, and z diiections) 
= 2.488 X 10*5 kN/m^ 

= 04 

= 2.7013 kN sVm'* 

= 20 X 10^ kN/m^ 

= 04 

= 1 0194 kN sVm^ 

= 75m 
= 2m 
= 10® kN 


= 100x125x50m 

= 30x10x15 (in x, y, and z directions) 
= 20 X 10® kN/m^ 

= 0.25 
= 16 kN /m® 

= 20 X 10“* kN/m^ 
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Poisson's Ratio Soil i/^ 

= 0.4 

Unit weight of Soil p. 

= lOkN/m^ 

Length of Pile / 

= 75m 

Diameter of Pile d. 

= 2m 

p 

= 10® kN 


w'' ^ Mz 



Figure 4 1 Triangular Pulse Excitation 


Displacement 
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Figure 4 4 Horizontal Displacement History of Pile Head Due to Triangular Pulse 
Load 


The displacement histories obtained with present model along with 
those of reported by Mamoon and Banerjee (1992) using time-domain 
BEM are presented in Figures 4,3, 4.4. It can be observed that the 
results obtained with present model compare well. 


4.5 Response of Single Piles 

Limited parametric studies have been conducted on the response of 
single piles. This has been guided by the fact that the similar behaviour 
of single pile and group of piles with a pile cap can not be expected. 
The compliance functions defined in Chapter 3 (Section 4) have been 
utilized for presentation of results by replacing the half-length L and 
half-width B of foundation block by the half the side d of the pile 
cross section. Compliance functions reported are computed using the 
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displacements at the pile head. 

4.5.1 Effect of Slenderness Ratio 

The slenderness ratio of the pile is defined as ratio of length of pile {/) 
and side (d) of square pile. For different l/d ratios pile head responses 
have been obtained in different modes of vibration. Vertical pile head 
1 espouse foi different l/d ratios for a given modular ratio Ep/Es are 
shown in Figures 4 5 and 4.6. It can be observed that, with increase 
in l/d ratio, compliance functions reduce, and reduction in comphance 
functions is more at low fiequency ratios. 



Figure 4.5: Vertical Compliance Function of Pile Head {Ep/Es — 10*) 
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Figure 4 6 Vertical Compliance Function of Pile Flead {Ep/Es = 10*') 


For horizontal mode of vibration, the compliance functions are shown 
in Figures 4.7 and 4.8. It can be seen that for low modular ratio 
[Ep/Ej = 10^) the l/d ratio does not have significant effect, whereas 
substantial reduction in compliance function with inciease in l/d ratio 
can be observed for modular ratios of 10'‘. 
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Frequency Ratio 


Figure 4.7. Horizontal Compliance Function of Pile Head {EpfEs = 10^) 



Figure 4 8; Horizontal Compliance Function of Pile Head [EpJEs = 10^) 

Since the piles are also modeled as continuum elements the rotational 
moment is simulated using a couple and the rotational displacement 
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is measured in terms slope of pile head. For rotational mode of 
vibration, the compliance functions are shown in Figures 4.9 and 4.10. 
The reduction in compliance functions with increase in slenderness 
ratio IS comparatively small for smaller modular ratio and high at 
higher modular ratios 



Figure 4 9- Rotational Compliance Function of Pile Head [EpjE, = 10^) 
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Figure 4 10' Rotational Compliance Function of Pile Head {E^fEs = 10'’) 


4.5.2 Effect of modular ratio 

The effect of modular ratio on the compliance functions for verti- 
cal/ horizontal and rotational modes are shown in Figures 4.11, 4.12 
and 4.13. It can be seen that modular ratio has significant effect at 
higher slenderness ratios. Vertical and rotational modes are strongly 
influenced by the slenderness ratio and modular ratios. 
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Figure 4.11 Vertical Compliance Functions of Pile Head for Different Modular 
Ratios 



Figure 4 12 Horizontal Compliance Functions of Pile Head for Different Modular 


Ratios 
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Figure 4 13: Rotational Compliance Functions of Pile Head for Different Modular 
Ratios 


4.6 Response of Pile Foundations 

In this section the response of pile foundation soil systems investi- 
gated are presented. The parameters that effect the response such 
as slenderness ratio of pileS/ number of piles ( i.e, spacing of piles) 
modular ratio of pile/pile cap material soil and effect of layered soils 
and thickness of pile cap, are considered for the study. The responses 
for different modes of vibration are presented using the compliance 
functions defined in Chapter 3 (Section 4). The geometry and the 
material properties of the system considered for illustration are given 
below. The pile foundation consists of pile cap of 1.0m thick and four 
pile at four corners of the cap. The arrangement of piles are shown 
in Figure 4,14 
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2B *\ 


(a) Pjifi roundiUnn wich 4 Pile* (b) Pile Foimdaiion wiA 8 Pifea 


Figure 4 14 Plan Showing the Arrangemciit of Piles 


Size of the model 
Number of elements 
Young's modulus of pile material Ep 
Poisson's ratio pile material i^p 
Unit weight of pile material pp 
Young's modulus of soil E, 

Poisson's ratio soil 
Unit weight of soil pj, 

Half-width of cap B 

Half-length of cap L 

Thickness of the cap t 

Depth of embeddment of the cap li 

Side of the square pile d 


= 100x50x50m 

= 30xl0xl5(in x, y, and z directions) 
= 47.03GPa 
= 0.33 
= 24.0kN/m3 
= 47.03MPa 
= 0.33 

= 17.0 kN/m^ 

= 2.0m 
= 4.0m 
= 1.0m 


= 1.0m 


= 10m 
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4.6.1 Effect of Slenderness Ratio 

The compliance funcUon of pile foundation with a rigid cap for dif- 
ferent modes are obtained with piles of varying slenderness ratios 
The vertical compliance of rigid cap with and without piles of dif- 
ferent slenderness latios for a particular modular ratios are shown m 
Figures 4.15 and 4.16. 



Figure 4 15. Vertical Compliance Functions of Pile Cap {EpjEg = 10^) 


I 

I 
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Figure 4 16 Vertical Compliance Functions of Pile Cap {EpfEj = 10'*) 


It can be observed that presence of piles reduce the vertical compliances 
by a factor of 4 to 5. The resonant frequencies reduce due to presence 
of piles and so is the sharpness at resonance. 

Horizontal compliance functions of the rigid cap of pile group are 
shown in Figures 4.17 and 4.18, The reduction of compliance functions 
due to piles is marginal (of the order of 50%), However the resonant 
frequencies reduce due to presence of piles. 
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Figure 4.17 Horizontal Compliance Functions of Pile Cap {Ep(E^ = 10^) 



Figure 4.18. Horizontal Compliance Functions of Pile Cap [EpJEs = 10“^) 


Rotational compliance functions of the pile cap for different l/d ratios 
arc shown in Figures 4,19 and 4.20. It can be observed that the 
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reduction in rotational compliance function due to the piesence of 
piles IS substantial. They reduce almost by a factor of 4 to 20. 



Figure 4 19 Rotational Compliance Functions of Pile Cap (/?,,//?, - 10'’) 



Figure 4 20. Rotational Compliance Functions of Pile Cap " i0'‘) 





4 6 Response of Pile Foundations 


145 


4.6.2 Effect of Modular Ratio 

The effect of niodulai ratio at different IJd ratios for different modes 
of vibration are shown m Figures 4.21, 4 22 and 4 23. It can be 
observed that Ihe increase in modular ratio does not significantly 
effect the vertical and horizontal compliances, whereas the rotational 
compliances decrease by a factor of 2 to 5. 



Figure 4.21: Effect of Modular Ratio on Vertical Compliance Functions of Pile 
Cap 



146 


PILE FOUNDATION SOIL SYSTEMS 



Figure 4 22; Effect of Modular Ratio on Horizontal Compliance Functions of Pile 
Cap 



Figure 4.23. Effect of Modular Ratio on Rotational Compliance Functions of Pile 
Cap 
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4.6.3 Effect of Layered Soil 

Compliance functions of pile foundations in two layered soil system 
have been investigated. The soil system considered has a top layer 
with varying thickness and the piles have a slenderness ratio of 20. 
The modular ratio M, of the two layers and the thickness ratio T; 
defined m the Chapter 4 (Section 4 5) has been used to piesent the 
results. The vertical, houzontal and rotational compliance functions 
of pile cap are shown for different thickness ratios in Figures 4 24, 
4.2.5 and 4.26. Presence of hard stratum underneath the soil reduces 
the compliance functions for all modes of vibration and incieases the 
resonant frequencies of the system. The reduction in compliance is 
highest for vertical and lowest for horizontal mode. Unlike block 
foundation, the pile foundation responses are influenced by soil of 
larger thickness. 



Figure 4 24: Effect of Layer Thickness on Vertical Compliance Functions of Pile 


Cap 
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Figure 4 25 Effect of Layer Thickness on Horizontal Compliance Functions of 
Pile Cap 



Figure 4.26: Effect of Layer Thickness on Rotational Compliance Functions of 
Pile Cap 
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4.6.4 liffect of Number of Piles 

'Jb invcslignte the effect of pile spacing/number of piles on the com- 
pliance functions, pile foundation with 8 piles is considered Here the 
piles are placed at middle of the sides in addition to the corners as 
shown in Figure 4.14 Vertical, horizontal and rotational comphance 
functions with 8 piles are shown in Figures 4 27, 4 28 and 4 29. The 
1 eduction in compliance due to reduced spacing is small for vertical 
and hoiizonlal modes whereas about 50% reduction is observed for 
rotational modes. 



Frequency Ratio 

Figure 4.27: Effect of Number of Piles on Vertical Compliance Functions of Pile 
Cap 
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I/d = 10, 4 Piles 
I'd = 10, 8 Piles 
I/d = 40, 8 Piles 
I/d = 40, 8 Piles 


Frequency Ratio 

Figure 4 28 Effect of Number of Piles on Horizontal Compliance Functions of 
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Figure 4.29; Effect of Number of Piles on Rotational Compliance Functions of 
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4.7 Response of Machine Pile Foundation Sys- 
tems With and Without Piles 

I'hc computnlional model developed has been utilized to analyze diesel 
gcncintoi foundation. A schematic diagram of the foundation system 
adopted is shown in Figure 4 32. The machine is transmitting lateral 
foice and rotational and torsional moments which are periodic but not 
hnimonic. The excitations are shown in Figures 4.30 and 4.31 The 
malenal properties and the geometry of the foundation system are 
given below. 


Size of the model 
Number of elements 
Young's modulus of pile material Ep 
Poisson's ratio pile material Up 
Unit weight of pile material Pp 
Young's modulus of soil 
Poisson's ratio soil 
Unit weight of soil p. 

Width of cap D 

Length of cap L 

Thickness of the cap t 

Depth of embeddment of the cap h 

Side of the square pile d 


30x20x30m 

20xl0xl5(in x, y, and z direchons) 

14.1GPa 

0.25 

24.0kN/m3 

20.0MPa 

0.25 

16.67 kN/m^ 

2.0m 

6.0m 

0.5m 


0.5m 

0.5m 



Force 
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Figure 4 30 Horizontal and Rotational Excitation 



Figure 4 31: Torsio 
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Figure 4 32. Machine Pile Foundation System 


Maximum displacements in each of the mode at the center of gravity of 
the pile cap are used to present the results. The computational model 
used is essentially the one described in Chapter 3. The numerical 
integration of the system of differential equations is carried out till 
the responses reach steady state conditions. The maximum responses 
(displacements) in steady state region are used to show the effect of 
pile cap thickness, effect of Iji ratio of the piles and the shear modulus 
of the soil. 

4.7.1 Effect of Shear Modulus 

The effect of shear modulus on the maximum displacement of the 
pile capA>lock in horizontal, rotational and torsional modes are shown 
in Figures 4.33, 4.34 and 4.35. The results shown are for a block 
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thickness of 0.5m and a slenderness ratio of 20 for the piles. It can 
be observed that increase in shear modulus beyond 6MPa does not a 
have significant effect on the responses. The presence of piles have 
greater effect in rotational and torsional modes than in horizontal 
mode. Shear modulus of the soil strongly influences the responses of 
block foundation compared to pile foundation. 



Figure 4.33: Effect of Shear Modulus on Horizontal Response 



Maximum Displacement (m) Maximum Displacement (m) 
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Figure 4 34. Effect of Shear Modulus on Rotational Response 



Figure 4.35' Effect of Shear Modulus on Torsional Response 
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4.7.2 Eifect of Block/Pile Cap Thickness 

The response of pile foundation system with shear modulus of lOMPa 
IS used to study the effects of slenderness ratio and pile cap/block 
thickness. The maximum displacements m horizontal, rotational and 
torsional modes for different values of thickness of pile cap/block, 
are shown in Figures 4.36, 4.37 and 4.38. Increase in pile cap/block 
thickness from 0.5m to 1.5m reduce the displacements marginally 
in horizontal mode where as the reduction is by a factor of 2 in 
torsional mode. Increase in thickness beyond 1.5m does not reduce 
the displacements in horizontal and rotational modes significantly. 
Increase m 1/d ratio does not significantly effect the responses. 



Figure 4 36 Effect of Pile Cap/Block Thickness on Horizontal Response 




Maximum DfspJacement (m) Maximum Displacement (m) 
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;iire 4.37: Effect of Pile Cap/Block Thickness on Rotational Response 



Figure 4.38- Effect of Pile Cap/Block Thickness on Torsional Response 
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4.8 Summary 

A brief summary of the observations made from the studies on the 
response of pile foundation system are presented in this section. 

1. The results obtained with present model agree well with those 
obtained with time-domam BEM. 

2. Slenderness ratio has a very significant effect on the vertical and 
rotational response of single piles. The reductions in compliances 
m these inodes are by a factor of 4-5 due to increase of Ijd ratio 
from 10 to 40. 

3. Modular ratio values beyond 10'‘ do not significantly effect the 
responses of single piles 

4. The presence of piles reduce the compliances by a factor of 2 to 
4 for vertical mode, and by less than 50% for horizontal mode 
and by a factor of 10 for rotational mode. 

5. Modular ratio does not have significant effect on the pile gioup 
response. 

6. The presence of piles reduce the resonant frequencies compared 
to the block foundation for all modes. 

7. Stiffer lower layer reduces the compliance functions of the pile 
group and increase the resonant frequency, Signihcant thickness 
of the top layer, effects the responses of the pile foundation in 
layered soil. 
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8 Increase in the number of piles reduces the compliances , but 
the 1 eduction is not pioporlional to the number of piles. Effect 
of highei number of piles is felt more at higher frequencies 

9. Inciease in shear modulus reduces the responses of pile founda- 
tion significantly at low modulii 

10. The pile cap/block thickness does not significantly effect the lateral 
and rotational responses, whereas it has significant effect on 
toisionni icsponse, 

The capabilities of the present computational model can be summarized 
as follows: 

1. The Piesent model can simulate semi-infinite half-space. 

2. The Present model can describe the dynamic behaviour of a three 
dimensional discretized model consisting of about 13,000 degrees 
of freedom by a small fraction of about 30 degrees of freedom. 

3. The present model can handle arbitrary excitations and the results 
obtained will be accurate, since continuous monitoring of local 
truncation error is carried out. 

4. The computational effort is limited, since the computadons are 
done with a few degrees of freedom. 

5. The present model can be applied to structures having regular 
geometry, non-homogeneous domains, non-convex domains 

6. An integrated approach wherein soil, structure, foundation can 
considered together without any need for substructuring is pre- 
sented in this model. 



Chapter 5 


ANALYSIS OF VIBRATION 
ISOLATION PROBLEMS 

5.1 Introduction 

Vibrations caused by machine foundations, traffic, blashng propagate 
through ground motion causing disturbances to adjacent structures 
and may affect the operation of sensitive equipment. Most of the 
energy affecting the adjacent structures is transmitted by Rayleigh 
waves which travel close to the surface of the ground (E E. Richart, 
Hall and Woods, 1970. Placing a wave barrier in the ground before 
the structure (receiver) to reduce the disturbances is called passive 
isolation. On the other hand placing a barrier near the source of 
vibration is called active isolation (E E. Richart, Hall and Woods, 1970) 
The effectiveness of the wave barriers depend upon the absorption or 
screening of the Rayleigh wave energy introduced by them. Generally 
open trenches, in-filled trenches, sheet pile walls, or rows of piles are 
provided to achieve vibration isolation. 
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In recent tunes vibration isolation is being studied using numerical 
and experimental research. Analytical or closed form solutions are 
extremely difficult, if not impossible, owing to the geometry and non 
homogeneity of the material involved An efficient numerical technique 
can provide an insight m to the vibrahon isolation phenomenon. 
Woods (1967) performed a series of field experiments on active isolation, 
where the trenches are put near the source, and passive isolation, 
where the trenches are placed at far held. Based on his experimental 
results, he presented a some guidelines for the dimensions of the 
open trench to achieve an amplitude reduction of 75% or more. In 
past couple of decades, number of researchers used various numerical 
techniques to study the vibration-screening problems. Aboudi (1973) 
used finite difference method (FDM) with special treatment of external 
boundaries to evaluate ground response of elastic half-space with thin 
barrier embedded. He concluded that better screening is observed 
m the case of shallow screen, by using low velocity and low density 
barriers. Haupt (1977) studied the vibration isolation of concrete walls. 
He used finite element method with an influence matrix boundary 
condition. He concluded that only cross section area influences the 
isolation effect but not the geometrical shape of the barrier. Segol, 
Lee and Abel (1978) developed a method in which vertical boundary 
forces were obtained m terms of modal matrix of eigen vectors and 
eigen values. They used finite element method for near field. They 
concluded that ratio of trench depth to shear wave length is an 
influencing factor of vibration isolation. Fuyuki and Matsumoto (1980) 
have investigated the Rayleigh wave scattering by rectangular trenches 
using finite difference method with improved absorbing boundary 
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conditions. Ihcy found that width of shallow open trenches could 
have sij;nlficant effect on vibration isolation May and Bolt (1982) have 
studied llie effectiveness of open tienches on vibration screening to 
hoiizontnlly propagating 1) SV and SH waves in a two-layered soil 
system. 

Recently^ boundary element methods have emerged as an alternahve 
numoiical scheme for solving wave propagation problems involving 
semi-infinile domains. Emad and Manolis (1985) have studied vibra- 
tion isolation by open trenches in frequency domain. Beskosy Dasgupta 
and Vnrdoulakis (1986) have studied vibration isolation of open and 
concrete trenches due to transient wave source in Laplace domain. 
Lcungy Beskos and Vardoulakis (1990) used thin layer Green's function 
developed by Kausel and Peek (1982) and Hull and Kausel (1984) to 
study the vibration isolation in layered soils. Dasgupta, Beskos and 
Vardoulakis (1988) and Banerjee, Ahmad and Chen (1988) applied 3Di- 
mensional BEM to vibration isolation by open trenches Beskos and 
Vardoulakis (1990) have presented results of earlier work on vibration 
screening of open and fiUed trenches. They used a frequency domain 
boundary element method for 3-D and plane strain idealization of 
waves generated by rigid vibrating foundations. They used full-space 
Green's functions for homogeneous soil and those developed by Kausel 
and Peek (1982) for layered soil. They showed that in-hlled trench 
performed better than open trench when 3— D analysis is adopted 
Chouw, Le and Schmid (1990) have studied vibration isolation of open 
trenches due to harmonically excited strip footings using boundary 
clement method. They have conducted a detailed parametric study 
on the effectiveness of the open trenches for vibration isolation due 
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to strip footing vibrating in different modes, Ahmad and Al-Hussnini 
(1991) have studied the vibration isolation by open and filled trenches 
using 2-Dimensional BEM using higher order elements in frequency 
domain They carried out parametric studies and developed some 
simplified models and suggested design guidelines. Al-Hussiani and 
Ahmad (1991) have studied vibration isolation by open trenches due 
to honzontal vibrating source. 

5.2 Present Model 

Most of the existing procedures for analysis of filled and open trenches 
use two dimensional modeling. Few of the three dimensional time 
domain BEMs use constant/linear elements. It has been pointed out 
(Von Estorff, Pais and Kausel, 1990) that the results obtained using 
time-domain methods do not obey the causality condition for non- 
convex domains like the present vertical trenches. As already pointed 
out the time-domain BEMs require huge computational resources. 
Most of the studies ignore the effect of trench width possibly due to 
the presumed fact that the ^3 of energy is transmitted by the surface 
waves. Recently Meek and Wolf (1993) have shown that the partillon 
of the energy between the body waves and surface waves change 
with frequency ratio, and at higher frequency ratios body waves carry 
higher energies. Another anomaly observed by Beskos, Dasgupta and 
Vardoulakis (1986) is that 2D modeling and 3D modeling resulted in 
opposite conclusions about the efficiency of open and filled trenches. 
There exist a need for systematic study of open and filled trenches 
using 3D models in time-domain. The effect of layering, the location. 
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wkHIi and depth of open and filled trench are some of the parameters 
that need detailed investigation in three dimensions The present 
study makes use of the model described in Chapter 2 to study these 
pioblems. The modeling consists of - 1. discretizing a ktgc domain 
using the 3D 8 noded elements and generating mass stiffness and 
damping matrices and load vector which describe equations of motion, 
2. Genciating load dependent Lanczos vectors and transforming the 
equations of motion to Lanczos space, and 3. Solving the transformed 
equations of motion is time domain by adaptive time stepping scheme 
and obtaining the responses in original space 


5.3 Presentation of The Results 

The problem definition of an active vibration isolation system in three 
dimensional idealization using an open and filled (barrier) trenches 
are shown in Figure 5 1. The geometry of the system (cross section 
across the footing) is shown in Figure 5,2. The dimensions of the 
Ircnch/barrier - Location /, depth d, width w and width of foundahon 
Wf of the Irencli/barrier are expressed in non-dimensional quantities 

as, 


L 

D ^Lr(t 

11 ^ ^LrW 

Wf =LrWf 


Since most of the energy 


(5.1) 

transmission takes place through Rayleigh 
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waves, and consequent screening of it lesults in vibration isolation, the 
depth, width and location of the trench/barrier are normalized with 
respect to Rayleigh wave length. The efficiency of the trench is usually 
measured in terms of amplitude reduction behind the trenclVbarner. 
A non-dimensional parameter A, (Amphtude reduction ratio) at any 
point, expressed as the ratio of amphtude of the ground surface with 
trench to that of without trencIVbarrier is used to express the vibration 
isolation efficiency Contours of amphtude reduction ratio /I, are used 
to present the results. The properties of the footing are used for 
material of the barrier also. 




Vt 


(5.2) 


where vt, and vt are the amplitudes with and without trench respec- 
tively. 


5.3.1 Comparison of The Model 

A machine foundation system surrounded by an open trench at shoit 
distance acting as an active isolation system studied by Dasgupta and 
Vardoulakis (1986) has been used for comparison. The geometry and 
the problem definition is shown in Figure 5.1. The material and 
geometric parameters of the system are given below. 


f) 3 Presentation of The Result 



Figure 5 1 Open 
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Size of the model 


12xl2xl2m 

Number of Elements 


18x10x18 

Young's Modulus of Footing Matenal Ep 


11315 MPa 

Poisson's Ratio Footing Material i/p 


0.25 

Unit Weight of Footing Material Pp 

— 

24.0 kN/m^ 

Young's Modulus of Soil 


330 MPa 

Poisson's Raho Soil i/j 


0.25 

Unit Weight of Soil pa 


17.0 

Width of the Trench w 


0.3m 

Depth of the Trench d 


2.5m 

Location of the Trench 1 


2 Om 


The contours of amplitude reduction ratio of the surface for the trench 
obtained with present model and the that of Dasgupta and Vardoulakis 
(1986) using linear elements and of Banerjee, Ahmad and Chen (1988) 
using higher order elements are shown in Figure 5.3. It can be seen 
that the difference between the two results is small behind the trench. 
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(c) Banerjee, Ahmad and Chen 
(1988) Results 

rigure 5 3: Comparison of Amplitude Reduction Ratio Contours 


5.3.2 Effect of Treiich/Barrier Depth 

Since effective screening of the wave energy leads to proper vibration 
isolation and Rayleigh waves carry most of it, the length of the path 
before the waves reach the target (area which is to be screened) 
determines the efficiency. An easy way to achieve such an objective 
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is to have trench with sufficient depth near the target. The non 
dirnensional trench/barrier depth D defined in (5.1) is used to study 
the effect. The contours of amplitude reduction ratio of the suiface 
for different values of D are shown in Figure 5.4 Successful vibration 
isolation is defined by F. E. Richarh Hall and Woods (1970) as the one 
which yield an amplitude reduction ratio of 0 25. It can be seen that 
the trench with a normalized depth of 0.27 does not screen effectively, 
yielding a reduction ratio of 0 7, where as the trench with D =- 0.82 
gives an amplitude reduction ratio of 0 2 over a large area. The open 
trench causes an amplification of the displacements in front of the 
trench to as high as 2. 

The screening effect of barrier has also been studied with different 
depths of barrier The amplitude reduction is shown in Figures 5.5 for 
different values of nondimensional depth parameter D, The amplitude 
reduction ratio achieved with a barrier is on the higher side compared 
to that of an open trench. The maximum reduction ratio achieved is 
of the order of 0 4 only with a barrier of depth D - 0.82. Flowever 
there is no amplification of the displacements in front of the bairier 
as in the case of the trench. 
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Figure 5,4. Effect of Trench Depth on Amplitude Reduction Ratio 


The screening effect of bairier has also been studied with different 
depths of barrier. The amplitude reduction is shown in Figures 5 5 for 
different values of nondimensional depth parameter D. The amplitude 
reduction ratio achieved with a barrier is on the higher side compared 
to that of an open trench. The maximum reduction ratio achieved is 
of the order of 0.4 only with a barrier of depth D = 0,82, however 
there is no amplification of the displacements in front of the barrier 
as in the case of the trench. 
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(a) D = 0 27 (b) D = 0 5 
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(c) D = 0.82 

Figure 5 5 Effect of Barrier Depth on Amplitude Reduction Ratio 

5.3,3 Effect of Trench/Barrier Width 

The effect of trench/barrier width on vibration isolation has been 
investigated. The amplitude reduction ratio contours for different 
nondimensional widths W for a depth ratio £) of 0 6 is shown in 
Figures 5.6 and 5,7. It can be seen that the amplitude reduction ratio 
IS not significantly influenced by the width of of the open trench. 
A trench width o( W ~ 0.6 gave higher amplitude reduction. In 
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case' of coiiciclc bariiei (Figure 5,7) the width of the barriei does not 
have significant effect on amplitude leduction ratio. The amplitude 
1 eduction (scicening) extends to a small area behind the barrier 
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Figure 5.6: Effect of Trench Width on Amplitude Reduction Ratio 
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Figure 57 Effect of Barrier Width on Amplitude Reduction Ratio 

5,3.4 Effect of Layered Soils 

Vibration isolation efficiency of an open trench in a two layered soil 
system has been investigated. The system considered in this study 
consists of a finite thick top layer on half-^space. The thickness of the 
top layer and the modulii ratio of top and bottom soils have been 
used to present the results. A non-dimensional top thickness and 
modular ratio M, are defined as 
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’J\ --T/A, 

M, =Ei,/Ei (5 3 ) 

where t is the thickness of top soil layer, E, arc young modulii of 
bottom and top soils respectively (Figure 5 2) 

The ampliUide reduction ratio contours of the surface for layered soils 
with different top layei thickness values as well as for half-space 
with trench depth Z) = 0 5 are shown along with the results for 
homogeneous half-space in Figures 5 8 and 5,9 

It can be seen that softei layer beneath stiff top layer, ; c f\f, < 1 0, 
(Figure 5.8} reduces the vibration isolation efficiency of the trench. 
Increasing the top layer thickness increases the efficiency. Whereas a 
stiff layer beneath a soft soil, ie. Mr >10, improves the vibration 
isolation efficiency of the bench (Figure 5 9) 



(a) 7; = 0.5 (b) T; = 1.0 

Figure 5,8 Effect of Top Layer Thickness on Amplitude Reduction Ratio for 


Mr = 0.5 
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(a) T, =05 



(b) Tr = 1.0 (c) Half-space 

Figure 5 9 Effect of Top Layer Thickness on Amplitude Reduction Ratio of an 
Open Trench with Mr = 2.0 


The amplitude reduccion contours for concrete trench in two layered 
system is shown in Figure 5,10. It can be observed that stiffer bottom 
layer improves the screening efficiency of the barrier. As the thickness 
of top layer reduces the efficiency improves, 
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Figure 5,10^ Effect of Top Layer Thickness on Amplitude Reduction Ratio of 
Barrier with = 2.0 


5.3.5 Effect of Frequency Ratio 

As inenlioned in the previous section the ratio of energy transmitted 
by the Rayleigh wave reduces with increasing frequency ratio, the 
effect of frequency ratio on the efficiency of the open trenches is 
investigated. Frequency ratio is changed by changing the dimensions 
of the foundation, keeping the excitation frequency same. The effect of 
frequency ratio for given frequency for different widths of the trench 
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IS investigated. Amplitude reduction ratio contours are shown for 
different frequency ratios in Figure 5.11. It can be observed that, for 
wider trenches {te.W = 0.18), there is small improvement in efficiency 
at higher frequency ratio. For practical range of frequency ratios this 
effect is not significant. 
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5.4 Summary 

Based on the limited parametric studies carried out the following 

observations can be made, 

1. The present computational scheme predicts a higher amplitude 
reduction ratios compared to those obtained with BEMs. 

2. In general open trenches produce better screening than concrete 
barriers. However no amplification of the displacements occurs 
in front of the trench due barriers, 

3. Ihench/Barrier depths have significant effect on the vibration iso- 
lation efficiency. An open trench with a depth Z) = 0 82 produces 
successful screening. 

4. Trench/Barrier width does not significantly effect the vibration 
isolation. Similar observations was made by Woods (1967). 

5. Layering has stronger influence on concrete barriers than open 
trenches, 

6. Frequency ratio does not significantly effect the vibration isolation. 



Chapter 6 


CONCLUSIONS AND SCOPE FOR 
FUTURE STUDY 


6.1 Introduction 

The results obtained in the present investigation for different class of 
problems have already been discussed at the end of the respective 
chapters. The advantages of the computational scheme and observa- 
tions made from the results obtained for various class of problems and 
performance of the scheme are summanzed in the following section 
Possible extensions of the study and scope are indicated in Section 6 3. 

6.2 Conclusions 

The computational scheme developed in this investigation is capable of 
analyzing wide class of time domain soil-structure interaction problems 
accurately. The advantages of the method, as brought out from the 
study, are summanzed below. 
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1. A large mesh coupled with Lanczos transformation scheme could 
represent semi-infinite soil mass, with out any need for artificial 
boundaries, for computational purposes, 

2. The present scheme can handle transient as well as periodic exci- 
tations containing arbitrary amplitudes and different frequencies. 

3 An adaptive strategy is incorporated for error control due to tem- 
poral discretization and coordinate transformation, thus avoiding 
problems of stability, convergence etc. . 

4. The error due to transformation decreases exponentially with 
increase in number of Lanczos vectors. The incremental com- 
putational cost per vector remains almost constant, and hence 
desired level of accuracy can be obtained. 

5. Even though the transformation matrices are load dependent, the 
same can be used in similar loading conditions. 

6. A few degrees of freedom (about 20-40) could simulate the be- 
haviour of a physical model represented by few thousands de- 
grees of freedom. 

7. The present model can effectively analyze both convex as well as 
non-convex domains. The solutions obey causality conditions. 

8. The computational scheme is local both in space and time. 

9. The present scheme demands moderate computational resources 
in terms of storage, but requires small amount of processor time. 
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10 In the piesent scheme a unified analytical framework is avail- 
able for soil-structure interaction problems considering the soil, 
structure and foundation as an integrated model. 

Compliance functions of embedded foundations are obtained, using 
three dimensional modeling, for different modes of vibration Effect 
of embeddment, aspect ratio, flexibility of the footing and layering 
of soil, linearly varying modulus and finite thickness of the strata 
on compliance of the footing have been presented The present 
investigation presented a methodology to study the realistic response 
of the foundation in both time and frequency domain Based on the 
numerical simulations earned out the following observations can be 
made: 

1. Reduction in compliance functions due to embeddment is more 
m horizontal and rotational modes than in vertical mode. Em- 
beddment does not alter the resonant frequency of the system in 
any significant manner. 

2. Increase in aspect ratio increases the compliance funchons in all 
modes. However the effect is small m rotational modes. The 
effect of aspect ratio is small at higher frequencies 

3. Compliance functions are strongly influenced by a small thickness 
of soil at the surface. 

4. The foundation flexibility effects vertical comphances strongly. 

5. Linearly increasing modulus increases the resonant frequencies 
and reduces the sharpness of the peaks at resonance. 
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6. Presence of bedrock at finite depth influences vertical compliances 
strongly compared to those of other modes 

Response of single piles and group of piles v/ith cap have been inves- 
tigated The parameters that effect the response such as slenderness 
ratio of pile, modular ratio of pile material, thickness of pile cap, num- 
ber of piles and layering of soil have been investigated. The results in 
terms of compliance functions have been presented and observations 
made from the results obtained are summaiized The present model 
has been applied to investigate a real life problem of foundation for 
machine generating periodic but non harmonic excitation. A brief 
parametric study of the machine foundation is carried out to bring 
out the parameters which effect the response Some of the important 
observations made from the study are summarized below: 

1 Slenderness ratio has significant effect on vertical compliances of 
the pile. The decrease can be seen to be of the order of 1/4 to 
1/5 due to increase in Ijd ratio from 10 to 40. 

2 Modular ratio values beyond lO** do not effect the response of 
single piles. 

3. The presence of piles reduce the compliances of the block/cap 
by a factor of 2 to 4 for vertical mode, and by less than 50% in 
horizontal mode and by a factor of 10 for rotational mode. 

4. The presence of piles reduce the resonant frequencies. 

5. The pile capA>lock thickness does not significantly effect the lateral 
and rotational responses, whereas torsional responses reduce with 
increase in thickness. 
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6. The reduction in the response of the machine foundation to the 
piesence of piles is more at lower shear modulus 

A brief study is carried out on vibration isolation performance of open 
trenches and concrete barriers. A machine foundation surrounded by a 
trencli/barrier has been used to inveshgate the efficiency of the isolation 
system. The effect of depth, width of the trench/barner, the frequency 
ratio of the system, layering of soil on the efficiency of the vibrahon 
isolation system have been investigated. The results in terms of 
amplitude reduction ratio contours have been presented. A summary 
of the observations made from the parametric study conducted are 
presented below, 

1. Open trenches provide better screening than concrete barriers. 

The open trenches amplify the displacements m front of the 
trench 

2, Open trench with depth ratio £> = 0 82 provide successful vibra- 
tion isolation, 

3, Efficiency of the isolation system increases with depth of trench/barrier. 

4, Width of trenclVbarrier does not significantly influence the effi- 
ciency, 

5. Frequency ratio does not significantly effect the vibration isolation 
efficiency over wide range of trench widths. 

6. Layering has stronger influence on concrete barriers than open 
trenches. 
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6.3 Scope for Future Study 

A wide class of linear time-domain soil stiucture interaction problems 
have been analyzed in the present investigation. The objective of 
the study is to provide a comprehensive affordable computational 
framework to analyze problems in time domain, A logical extensions 
of the present study could be in the following directions. 

1. Incorporation of non-linear material models. 

2. Separate material model/elements to account for soil-structure 
interface. 

3. Effect of partial uplift due to impact, strong ground motion etc . 

4. Kinematic interaction of adjacent structures. 

5. More detailed investigation of active and passive isolation prob- 
lems. 

6. Efficiency of isolation systems against earthquakes (transient ex- 
citation). 
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